From 802a4416b46bb8f41830e0e85e1b4d658760742f Mon Sep 17 00:00:00 2001 From: Alexander Schmidt Date: Mon, 10 Feb 2014 17:21:56 +0100 Subject: [PATCH] Migrated sim algorithm. --- .gitmodules | 3 + latlib | 1 + u1casc-ordinary/.gitignore | 6 ++ u1casc-ordinary/CMakeLists.txt | 20 ++++ u1casc-ordinary/FindGSL.cmake | 135 ++++++++++++++++++++++++++ u1casc-ordinary/latlib | 1 + u1casc-ordinary/sim.hpp | 141 ++++++++++++++++++++++++++++ u1casc-ordinary/u1casc-ordinary.cpp | 77 +++++++++++++++ 8 files changed, 384 insertions(+) create mode 100644 .gitmodules create mode 160000 latlib create mode 100644 u1casc-ordinary/.gitignore create mode 100644 u1casc-ordinary/CMakeLists.txt create mode 100644 u1casc-ordinary/FindGSL.cmake create mode 120000 u1casc-ordinary/latlib create mode 100644 u1casc-ordinary/sim.hpp create mode 100644 u1casc-ordinary/u1casc-ordinary.cpp diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..6386119 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "latlib"] + path = latlib + url = treefish.org:~/public_git/phys/latlib.git diff --git a/latlib b/latlib new file mode 160000 index 0000000..c08c4b7 --- /dev/null +++ b/latlib @@ -0,0 +1 @@ +Subproject commit c08c4b78cd61f072d566aefaa872c5348ab129a0 diff --git a/u1casc-ordinary/.gitignore b/u1casc-ordinary/.gitignore new file mode 100644 index 0000000..332a25a --- /dev/null +++ b/u1casc-ordinary/.gitignore @@ -0,0 +1,6 @@ +CMakeFiles +*~ +Makefile +cmake_install.cmake +CMakeCache.txt +u1casc-ordinary diff --git a/u1casc-ordinary/CMakeLists.txt b/u1casc-ordinary/CMakeLists.txt new file mode 100644 index 0000000..40225bd --- /dev/null +++ b/u1casc-ordinary/CMakeLists.txt @@ -0,0 +1,20 @@ +cmake_minimum_required(VERSION 2.8) + +project(u1casc) + +# #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-ordinary u1casc-ordinary.cpp) +target_link_libraries(u1casc-ordinary o815 ${GSL_LIBRARIES} lat_neigh) diff --git a/u1casc-ordinary/FindGSL.cmake b/u1casc-ordinary/FindGSL.cmake new file mode 100644 index 0000000..486be73 --- /dev/null +++ b/u1casc-ordinary/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-ordinary/latlib b/u1casc-ordinary/latlib new file mode 120000 index 0000000..7076779 --- /dev/null +++ b/u1casc-ordinary/latlib @@ -0,0 +1 @@ +../latlib \ No newline at end of file diff --git a/u1casc-ordinary/sim.hpp b/u1casc-ordinary/sim.hpp new file mode 100644 index 0000000..af989c9 --- /dev/null +++ b/u1casc-ordinary/sim.hpp @@ -0,0 +1,141 @@ +#ifndef SIM_HPP +#define SIM_HPP + +#include +#include +#include + +#include "latlib/neigh.h" + +#define EPSILONU 1 +#define EPSILONPHI 0.5 + +class sim : public o815::sim { +public: + sim(o815 *_O815); + unsigned int lsize4; + neigh *nb; + complex *U, *phi; + double kappa[2], lambda[2], beta; + +private: + void _makeSweep(); + void _newParas(); + gsl_rng* rangsl; + double rhoPhi(const int& iphi, const int& x0, const complex& candPhi); + double rhoU(const int& x0, const int& nu0, const complex& candU); + int updatePhi(const int& iphi, const int& x0); + int updateU(const int& x0, const int& nu0); +}; + +sim::sim(o815 *_O815) : o815::sim( _O815, + sizeof(complex)* + _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2+4) ) { + + 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]); + + phi = (complex*)confMem; + U = (complex*)(confMem + sizeof(complex)*lsize4*2); + + rangsl = gsl_rng_alloc(gsl_rng_ranlxs0); + gsl_rng_set(rangsl, time(NULL)); +} + +void sim::_makeSweep() { + for( int ix=0; ixparaQ)["kappaone"]; + kappa[1] = (*O815->paraQ)["kappatwo"]; + lambda[0] = (*O815->paraQ)["lambdaone"]; + lambda[1] = (*O815->paraQ)["lambdatwo"]; + beta = (*O815->paraQ)["beta"]; + + for(int ix=0; ix candU = U[x0*4+nu0] * polar(1.0, 2*EPSILONU*( 0.5 - gsl_rng_uniform(rangsl) )); + + if ( gsl_rng_uniform(rangsl) <= rhoU(x0, nu0, candU) ) { + U[x0*4 + nu0] = candU; + return 1; + } + + return 0; +} + +int sim::updatePhi(const int& iphi, const int& x0) +{ + complex candPhi = phi[ iphi*lsize4 + x0 ] + + complex(2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ), + 2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) )); + + if ( gsl_rng_uniform(rangsl) <= rhoPhi(iphi, x0, candPhi) ) { + phi[ iphi*lsize4 + x0 ] = candPhi; + return 1; + } + + return 0; +} + +double sim::rhoPhi(const int& iphi, const int& x0, const complex& candPhi) +{ + double deltaS=0; + + for( int mu=0; mu<4; mu++) { + if( iphi == 0 ) { + deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] ); + deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] ); + deltaS -= 2 * real( conj(candPhi) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] ); + deltaS -= 2 * real( conj(candPhi) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] ); + } + else if( iphi == 1 ) { + deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] ); + deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] ); + deltaS -= 2 * real( conj(candPhi) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] ); + deltaS -= 2 * real( conj(candPhi) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] ); + } + } + + deltaS -= kappa[iphi] * norm(phi[ iphi*lsize4 + x0 ]); + deltaS += kappa[iphi] * norm(candPhi); + + deltaS -= lambda[iphi] * pow(norm(phi[ iphi*lsize4 + x0 ]),2); + deltaS += lambda[iphi] * pow(norm(candPhi),2); + + return exp(-deltaS); +} + +double sim::rhoU(const int& x0, const int& nu0, const complex& candU) +{ + double deltaS=0; + + for( int nu=0; nu<4; nu++ ) { + if( nu == nu0 ) continue; + deltaS += beta * real( U[x0*4+nu0] * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) ); + deltaS += beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(U[ x0*4 + nu0 ]) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) ); + deltaS -= beta * real( candU * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) ); + deltaS -= beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(candU) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) ); + } + + deltaS += 2 * real( conj(phi[ 0*lsize4 + x0 ]) * U[ x0*4 + nu0 ] * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ] ); + deltaS -= 2 * real( conj(phi[ 0*lsize4 + x0 ]) * candU * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ] ); + + deltaS += 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(U[ x0*4 + nu0 ]) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ] ); + deltaS -= 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(candU) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ] ); + + return exp(-deltaS); +} + +#endif diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp new file mode 100644 index 0000000..f6f05fd --- /dev/null +++ b/u1casc-ordinary/u1casc-ordinary.cpp @@ -0,0 +1,77 @@ +#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" + +o815::comoption specOps[] = { + { "kappaone", required_argument, NULL, 'k', "set inverse mass kappa_1", "min:max:inc" }, + { "kappatwo", required_argument, NULL, 'm', "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" }, + { "", 0, NULL, 0, "", "" } +}; + +void parseSpecOps() +{ + for (int isopt = 0; isopt < O815->parsedSpecOps.size(); isopt++) + switch(O815->parsedSpecOps[isopt].first) { + case 'm': + O815->paraQ->addRange("mass", O815->parsedSpecOps[isopt].second); + break; + } +} + +void helpHeader() +{ + cout << "Usage: ./u1casc-ordinary [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, "phi2_hist") == 0 ) { + *O815->out->log << "MASTER: registered observable: phi2_hist" << endl << flush; + O815->observables.push_back(new obs_phi2_hist(O815)); + } + else if ( strcmp(*lonit, "phip2_hist") == 0 ) { + *O815->out->log << "MASTER: registered observable: phip2_hist" << endl << flush; + O815->observables.push_back(new obs_phip2_hist(O815)); + } + */ + } +} + +int main (int argc, char *argv[]) +{ + O815 = new o815(argc, argv, "u1casc-ordinary", specOps, &helpHeader); + + O815->addPara("beta", 1); + O815->addPara("lambdaone", 1); + O815->addPara("lambdatwo", 1); + O815->addPara("kappaone", 1); + O815->addPara("kappatwo", 1); + + parseSpecOps(); + O815->postParaInit(); + O815->Sim = new sim(O815); + parseLonelyArgs(); + O815->mainLoop(); + + delete O815; + return 0; +} -- 2.39.5