Performing valid state test. Added prestore.

diff --git a/.gitignore b/.gitignore
index 7893817f65abc1b7d9b659d42f02a47cd3b497f1..f898742879d543f40804765cf49304e0cd83f6a9 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,4 @@ CMakeFiles
 cmake_install.cmake
 Makefile
 clebschgordan
+su2clebschgordan

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 70efbe2911a331d891a37df9a4357089a7a22745..8548512d3254e2235cd58a5c376ccb6fb31b6d9a 100644 (file)
--- 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 49a80478de07aec77834719deb2c56107d1a09c8..11d506e5b1667863d24a21c147b84deb5a7c1dd7 100644 (file)
--- 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 689c0a1424fb1e7b94902b4749270081a4d5b563..87e01bf75020b7d2ffc5414e53ce3d196edc2d73 100644 (file)
--- 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 (file)
index 0000000..968c9de
--- /dev/null
+++ b/su2clebschgordan.cpp
@@ -0,0 +1,44 @@
+#include <iostream>
+
+#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;
+}