clebsch.h

   1 #ifndef CLEBSCH_H
   2 #define CLEBSCH_H
   3
   4 //#define NDEBUG
   5
   6 #include <algorithm>
   7 #include <cassert>
   8 #include <cmath>
   9 #include <cstdio>
  10 #include <cstdlib>
  11 #include <cstring>
  12 #include <functional>
  13 #include <fstream>
  14 #include <iostream>
  15 #include <map>
  16 #include <numeric>
  17 #include <vector>
  18
  19 // Declaration of LAPACK subroutines
  20 // Make sure the data types match your version of LAPACK
  21
  22 extern "C" void dgesvd_(char const* JOBU,
  23                         char const* JOBVT,
  24                         int const* M,
  25                         int const* N,
  26                         double* A,
  27                         int const* LDA,
  28                         double* S,
  29                         double* U,
  30                         int const* LDU,
  31                         double* VT,
  32                         int const* LDVT,
  33                         double* WORK,
  34                         int const* LWORK,
  35                         int *INFO);
  36
  37 extern "C" void dgels_(char const* TRANS,
  38                        int const* M,
  39                        int const* N,
  40                        int const* NRHS,
  41                        double* A,
  42                        int const* LDA,
  43                        double* B,
  44                        int const* LDB,
  45                        double* WORK,
  46                        int const* LWORK,
  47                        int *INFO);
  48
  49 namespace clebsch {
  50     const double EPS = 1e-12;
  51
  52     // binomial coefficients
  53     class binomial_t {
  54         std::vector<int> cache;
  55         int N;
  56
  57     public:
  58         int operator()(int n, int k);
  59     };
  60
  61     // Eq. (19) and (25)
  62     class weight {
  63         std::vector<int> elem;
  64
  65     public:
  66         // the N in "SU(N)"
  67         const int N;
  68
  69         // create a non-initialized weight
  70         weight(int N);
  71
  72         // create irrep weight of given index
  73         // Eq. (C2)
  74         weight(int N, int index);
  75
  76         // assign from another instance
  77         clebsch::weight &operator=(const clebsch::weight &w);
  78
  79         // access elements of this weight (k = 1, ..., N)
  80         int &operator()(int k);
  81         const int &operator()(int k) const;
  82
  83         // compare weights
  84         // Eq. (C1)
  85         bool operator<(const weight &w) const;
  86         bool operator==(const weight &w) const;
  87
  88         // element-wise sum of weights
  89         clebsch::weight operator+(const weight &w) const;
  90
  91         // returns the index of this irrep weight (index = 0, 1, ...)
  92         // Eq. (C2)
  93         int index() const;
  94
  95         // returns the dimension of this irrep weight
  96         // Eq. (22)
  97         long long dimension() const;
  98     };
  99
 100     // Eq. (20)
 101     class pattern {
 102         std::vector<int> elem;
 103
 104     public:
 105         // the N in "SU(N)"
 106         const int N;
 107
 108         // copy constructor
 109         pattern(const pattern &pat);
 110
 111         // create pattern of given index from irrep weight
 112         // Eq. (C7)
 113         pattern(const weight &irrep, int index = 0);
 114
 115         // access elements of this pattern (l = 1, ..., N; k = 1, ..., l)
 116         int &operator()(int k, int l);
 117         const int &operator()(int k, int l) const;
 118
 119         // find succeeding/preceding pattern, return false if not possible
 120         // Eq. (C9)
 121         bool operator++();
 122         bool operator--();
 123
 124         // returns the pattern index (index = 0, ..., dimension - 1)
 125         // Eq. (C7)
 126         int index() const;
 127
 128         // returns the pattern weight
 129         // Eq. (25)
 130         clebsch::weight get_weight() const;
 131
 132         // returns matrix element of lowering operator J^(l)_-
 133         // between this pattern minus M^(k,l) and this pattern
 134         // (l = 1, ..., N; k = 1, ..., l)
 135         // Eq. (28)
 136         double lowering_coeff(int k, int l) const;
 137
 138         // returns matrix element of raising operator J^(l)_+
 139         // between this pattern plus M^(k,l) and this pattern
 140         // (l = 1, ..., N; k = 1, ..., l)
 141         // Eq. (29)
 142         double raising_coeff(int k, int l) const;
 143     };
 144
 145     class decomposition {
 146         std::vector<clebsch::weight> weights;
 147         std::vector<int> multiplicities;
 148
 149     public:
 150         // the N in "SU(N)"
 151         const int N;
 152
 153         // save given irreps for later use
 154         const weight factor1, factor2;
 155
 156         // construct the decomposition of factor1 times factor2 into irreps
 157         // Eq. (31)
 158         decomposition(const weight &factor1, const weight &factor2);
 159
 160         // return the number of occurring irreps
 161         int size() const;
 162
 163         // access the occurring irreps
 164         // j = 0, ..., size() - 1
 165         const clebsch::weight &operator()(int j) const;
 166
 167         // return the outer multiplicity of irrep in this decomposition
 168         int multiplicity(const weight &irrep) const;
 169     };
 170
 171     class index_adapter {
 172         std::vector<int> indices;
 173         std::vector<int> multiplicities;
 174
 175     public:
 176         // the N in "SU(N)"
 177         const int N;
 178
 179         // save given irreps for later use
 180         const int factor1, factor2;
 181
 182         // construct this index_adapter from a given decomposition
 183         index_adapter(const clebsch::decomposition &decomp);
 184
 185         // return the number of occurring irreps
 186         int size() const;
 187
 188         // access the occurring irreps
 189         int operator()(int j) const;
 190
 191         // return the outer multiplicity of irrep in this decomposition
 192         int multiplicity(int irrep) const;
 193     };
 194
 195     class coefficients {
 196         std::map<std::vector<int>, double> clzx;
 197
 198         // access Clebsch-Gordan coefficients in convenient manner
 199         void set(int factor1_state,
 200                  int factor2_state,
 201                  int multiplicity_index,
 202                  int irrep_state,
 203                  double value);
 204
 205         // internal functions, doing most of the work
 206         void highest_weight_normal_form(); // Eq. (37)
 207         void compute_highest_weight_coeffs(); // Eq. (36)
 208         void compute_lower_weight_coeffs(int multip_index, int state, std::vector<char> &done); // Eq. (40)
 209
 210     public:
 211         // the N in "SU(N)"
 212         const int N;
 213
 214         // save irreps and their dimensions for later use
 215         const weight factor1, factor2, irrep;
 216         const int factor1_dimension, factor2_dimension, irrep_dimension;
 217
 218         // outer multiplicity of irrep in this decomposition
 219         const int multiplicity;
 220
 221         // construct all Clebsch-Gordan coefficients of this decomposition
 222         coefficients(const weight &irrep, const weight &factor1, const weight &factor2);
 223
 224         // access Clebsch-Gordan coefficients (read-only)
 225         // multiplicity_index = 0, ..., multiplicity - 1
 226         // factor1_state = 0, ..., factor1_dimension - 1
 227         // factor2_state = 0, ..., factor2_dimension - 1
 228         // irrep_state = 0, ..., irrep_dimension
 229         double operator()(int factor1_state,
 230                           int factor2_state,
 231                           int multiplicity_index,
 232                           int irrep_state) const;
 233     };
 234 };
 235
 236 #endif