From: Alexander Schmidt Date: Wed, 27 Aug 2014 22:40:56 +0000 (+0200) Subject: Performing valid state test. Added prestore. X-Git-Url: http://git.treefish.org/~alex/phys/su2clebsch.git/commitdiff_plain/a02d886eef0f2fb5c434d094ff34064aa62eed2e Performing valid state test. Added prestore. --- diff --git a/.gitignore b/.gitignore index 7893817..f898742 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ CMakeFiles cmake_install.cmake Makefile clebschgordan +su2clebschgordan diff --git a/CMakeLists.txt b/CMakeLists.txt index 70efbe2..8548512 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,3 +18,6 @@ target_link_libraries(clebschgordan clebsch) add_library(su2clebsch su2clebsch.cpp) target_link_libraries(su2clebsch clebsch) + +add_executable(su2clebschgordan su2clebschgordan.cpp) +target_link_libraries(su2clebschgordan su2clebsch) diff --git a/su2clebsch.cpp b/su2clebsch.cpp index 49a8047..11d506e 100644 --- a/su2clebsch.cpp +++ b/su2clebsch.cpp @@ -7,35 +7,83 @@ using namespace std; namespace su2clebsch { + int ipow(int base, int exp) + { + int result = 1; + while (exp) + { + if (exp & 1) + result *= base; + exp >>= 1; + base *= base; + } + + return result; + } + + prestore::prestore(unsigned int _max2j) + { + max2j = _max2j; + cgcdata = new double[ipow(max2j+1,6)]; + + for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++) + for (int factor1_mp=0; factor1_mp<=max2j; factor1_mp++) + for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++) + for (int factor2_mp=0; factor2_mp<=max2j; factor2_mp++) + for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++) + for (int irrep_mp=0; irrep_mp<=max2j; irrep_mp++) + { + cgcdata[ factor1_2j*ipow(max2j+1,5) + factor1_mp*ipow(max2j+1,4) + + factor2_2j*ipow(max2j+1,3) + factor2_mp*ipow(max2j+1,2) + + irrep_2j*(max2j+1) + irrep_mp ] + = + cgc( factor1_2j, 2*factor1_mp-max2j, factor2_2j, 2*factor2_mp-max2j, + irrep_2j, 2*irrep_mp-max2j ); + } + } + + bool stateValid(int state_2j, int state_2m) + { + if ( abs(state_2m) > state_2j || (state_2j + state_2m) % 2 != 0 ) + return false; + else + return true; + } + + double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, + int irrep_2j, int irrep_2m, const prestore& Prestore) + { + if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) ) + return 0.0; + + if ( factor1_2j > Prestore.max2j || factor2_2j > Prestore.max2j || irrep_2j > Prestore.max2j ) + return cgc(factor1_2j, factor1_2m, factor2_2j, factor2_2m, irrep_2j, irrep_2m); + + return Prestore.cgcdata[ factor1_2j*ipow(Prestore.max2j+1,5) + + (factor1_2m+Prestore.max2j)/2 *ipow(Prestore.max2j+1,4) + + factor2_2j*ipow(Prestore.max2j+1,3) + + (factor2_2m+Prestore.max2j)/2 *ipow(Prestore.max2j+1,2) + + irrep_2j*(Prestore.max2j+1) + + (irrep_2m+Prestore.max2j)/2 ]; + } double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, int irrep_2j, int irrep_2m) { + if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) ) + return 0.0; + const int factor1_n = (factor1_2j + factor1_2m) / 2; const int factor2_n = (factor2_2j + factor2_2m) / 2; const int irrep_n = (irrep_2j + irrep_2m) / 2; clebsch::weight factor1_S(2); factor1_S(1) = factor1_2j; - if ( ! (0 <= factor1_n && factor1_n < factor1_S.dimension()) ) { - //std::cerr << "clebsch warning: 0 <= factor1_state && factor1_state < factor1_dimension" << std::endl; - return 0.0; - } - clebsch::weight factor2_S(2); factor2_S(1) = factor2_2j; - if ( ! (0 <= factor2_n && factor2_n < factor2_S.dimension()) ) { - //std::cerr << "clebsch warning: 0 <= factor2_state && factor2_state < factor2_dimension" << std::endl; - return 0.0; - } - clebsch::weight irrep_S(2); irrep_S(1) = irrep_2j; - if ( ! (0 <= irrep_n && irrep_n < irrep_S.dimension()) ) { - //std::cerr << "clebsch warning: 0 <= irrep_state && irrep_state < irrep_dimension" << std::endl; - return 0.0; - } - + clebsch::decomposition decomp(factor1_S, factor2_S); bool b = true; @@ -46,8 +94,6 @@ namespace su2clebsch { } } if (b) { - cerr << "WARNING: 2j=" << irrep_2j - << " does not occour in decomposition!" << endl; return 0.0; } diff --git a/su2clebsch.hpp b/su2clebsch.hpp index 689c0a1..87e01bf 100644 --- a/su2clebsch.hpp +++ b/su2clebsch.hpp @@ -3,9 +3,22 @@ namespace su2clebsch { + struct prestore { + private: + unsigned int max2j; + double *cgcdata; + public: + prestore(unsigned int _max2j); + + friend double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, + int irrep_2j, int irrep_2m, const prestore& Prestore); + }; + double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, int irrep_2j, int irrep_2m); + double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, + int irrep_2j, int irrep_2m, const prestore& Prestore); }; #endif diff --git a/su2clebschgordan.cpp b/su2clebschgordan.cpp new file mode 100644 index 0000000..968c9de --- /dev/null +++ b/su2clebschgordan.cpp @@ -0,0 +1,44 @@ +#include + +#include "su2clebsch.hpp" + +using namespace std; + +int main () +{ + int max2j=10; + + cout << "Prestoring..." << flush; + su2clebsch::prestore su2cstore(max2j); + cout << "DONE" << endl; + + cout << "Checking..." << flush; + for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++) + for (int factor1_mp=0; factor1_mp<=max2j; factor1_mp++) + for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++) + for (int factor2_mp=0; factor2_mp<=max2j; factor2_mp++) + for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++) + for (int irrep_mp=0; irrep_mp<=max2j; irrep_mp++) + if ( su2clebsch::cgc( factor1_2j, 2*factor1_mp-max2j, factor2_2j, 2*factor2_mp-max2j, + irrep_2j, 2*irrep_mp-max2j, su2cstore ) + - + su2clebsch::cgc( factor1_2j, 2*factor1_mp-max2j, factor2_2j, 2*factor2_mp-max2j, + irrep_2j, 2*irrep_mp-max2j ) + != 0 ) { + cout << "CGC check failed!" << endl; + exit(1); + } + cout << "OK" << endl; + + cout << "Speedtest..." << flush; + for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++) + for (int factor1_mp=0; factor1_mp<=max2j; factor1_mp++) + for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++) + for (int factor2_mp=0; factor2_mp<=max2j; factor2_mp++) + for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++) + for (int irrep_mp=0; irrep_mp<=max2j; irrep_mp++) + double test = su2clebsch::cgc( factor1_2j, 2*factor1_mp-max2j, factor2_2j, 2*factor2_mp-max2j, + irrep_2j, 2*irrep_mp-max2j, su2cstore ); + + cout << "DONE" << endl; +}