From: Alexander Schmidt Date: Tue, 26 Aug 2014 22:33:13 +0000 (+0200) Subject: initial commit. X-Git-Url: http://git.treefish.org/~alex/phys/su2clebsch.git/commitdiff_plain/e2c2ab8026276b7ebe343b1ebe86f3248547940a initial commit. --- e2c2ab8026276b7ebe343b1ebe86f3248547940a diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7893817 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +CMakeCache.txt +CMakeFiles +*~ +*.a +cmake_install.cmake +Makefile +clebschgordan diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..70efbe2 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,20 @@ +project(SU2CLEBSCH) + +cmake_minimum_required(VERSION 3.0) + +set(CMAKE_BUILD_TYPE Release) + +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + +#set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}) + +find_package(LAPACK REQUIRED) + +add_library(clebsch clebsch.cpp) +target_link_libraries(clebsch lapack) + +add_executable(clebschgordan clebschgordan.cpp) +target_link_libraries(clebschgordan clebsch) + +add_library(su2clebsch su2clebsch.cpp) +target_link_libraries(su2clebsch clebsch) diff --git a/clebsch.cpp b/clebsch.cpp new file mode 100644 index 0000000..f5f852d --- /dev/null +++ b/clebsch.cpp @@ -0,0 +1,816 @@ +#include "clebsch.h" + +// implementation of "binomial_t" starts here + +clebsch::binomial_t binomial; + +int clebsch::binomial_t::operator()(int n, int k) { + if (N <= n) { + for (cache.resize((n + 1) * (n + 2) / 2); N <= n; ++N) { + cache[N * (N + 1) / 2] = cache[N * (N + 1) / 2 + N] = 1; + for (int k = 1; k < N; ++k) { + cache[N * (N + 1) / 2 + k] = cache[(N - 1) * N / 2 + k] + + cache[(N - 1) * N / 2 + k - 1]; + } + } + } + + return cache[n * (n + 1) / 2 + k]; +} + +// implementation of "weight" starts here + +clebsch::weight::weight(int N) : elem(N), N(N) {} + +clebsch::weight::weight(int N, int index) : elem(N, 0), N(N) { + for (int i = 0; index > 0 && i < N; ++i) { + for (int j = 1; binomial(N - i - 1 + j, N - i - 1) <= index; j <<= 1) { + elem[i] = j; + } + + for (int j = elem[i] >> 1; j > 0; j >>= 1) { + if (binomial(N - i - 1 + (elem[i] | j), N - i - 1) <= index) { + elem[i] |= j; + } + } + + index -= binomial(N - i - 1 + elem[i]++, N - i - 1); + } +} + +clebsch::weight &clebsch::weight::operator=(const clebsch::weight &w) { + int &n = const_cast(N); + elem = w.elem; + n = w.N; + return *this; +} + +int &clebsch::weight::operator()(int k) { + assert(1 <= k && k <= N); + return elem[k - 1]; +} + +const int &clebsch::weight::operator()(int k) const { + assert(1 <= k && k <= N); + return elem[k - 1]; +} + +bool clebsch::weight::operator<(const weight &w) const { + assert(w.N == N); + for (int i = 0; i < N; ++i) { + if (elem[i] - elem[N - 1] != w.elem[i] - w.elem[N - 1]) { + return elem[i] - elem[N - 1] < w.elem[i] - w.elem[N - 1]; + } + } + return false; +} + +bool clebsch::weight::operator==(const weight &w) const { + assert(w.N == N); + + for (int i = 1; i < N; ++i) { + if (w.elem[i] - w.elem[i - 1] != elem[i] - elem[i - 1]) { + return false; + } + } + + return true; +} + +clebsch::weight clebsch::weight::operator+(const weight &w) const { + weight result(N); + + transform(elem.begin(), elem.end(), w.elem.begin(), result.elem.begin(), std::plus()); + + return result; +} + +int clebsch::weight::index() const { + int result = 0; + + for (int i = 0; elem[i] > elem[N - 1]; ++i) { + result += binomial(N - i - 1 + elem[i] - elem[N - 1] - 1, N - i - 1); + } + + return result; +} + +long long clebsch::weight::dimension() const { + long long numerator = 1, denominator = 1; + + for (int i = 1; i < N; ++i) { + for (int j = 0; i + j < N; ++j) { + numerator *= elem[j] - elem[i + j] + i; + denominator *= i; + } + } + + return numerator / denominator; +} + +// implementation of "pattern" starts here + +clebsch::pattern::pattern(const pattern &p) : elem(p.elem), N(p.N) {} + +clebsch::pattern::pattern(const weight &irrep, int index) : + elem((irrep.N * (irrep.N + 1)) / 2), N(irrep.N) { + for (int i = 1; i <= N; ++i) { + (*this)(i, N) = irrep(i); + } + + for (int l = N - 1; l >= 1; --l) { + for (int k = 1; k <= l; ++k) { + (*this)(k, l) = (*this)(k + 1, l + 1); + } + } + + while (index-- > 0) { + bool b = ++(*this); + + assert(b); + } +} + +int &clebsch::pattern::operator()(int k, int l) { + return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1]; +} + +const int &clebsch::pattern::operator()(int k, int l) const { + return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1]; +} + +bool clebsch::pattern::operator++() { + int k = 1, l = 1; + + while (l < N && (*this)(k, l) == (*this)(k, l + 1)) { + if (--k == 0) { + k = ++l; + } + } + + if (l == N) { + return false; + } + + ++(*this)(k, l); + + while (k != 1 || l != 1) { + if (++k > l) { + k = 1; + --l; + } + + (*this)(k, l) = (*this)(k + 1, l + 1); + } + + return true; +} + +bool clebsch::pattern::operator--() { + int k = 1, l = 1; + + while (l < N && (*this)(k, l) == (*this)(k + 1, l + 1)) { + if (--k == 0) { + k = ++l; + } + } + + if (l == N) { + return false; + } + + --(*this)(k, l); + + while (k != 1 || l != 1) { + if (++k > l) { + k = 1; + --l; + } + + (*this)(k, l) = (*this)(k, l + 1); + } + + return true; +} + +int clebsch::pattern::index() const { + int result = 0; + + for (pattern p(*this); --p; ++result) {} + + return result; +} + +clebsch::weight clebsch::pattern::get_weight() const { + clebsch::weight result(N); + + for (int prev = 0, l = 1; l <= N; ++l) { + int now = 0; + + for (int k = 1; k <= l; ++k) { + now += (*this)(k, l); + } + + result(l) = now - prev; + prev = now; + } + + return result; +} + +double clebsch::pattern::lowering_coeff(int k, int l) const { + double result = 1.0; + + for (int i = 1; i <= l + 1; ++i) { + result *= (*this)(i, l + 1) - (*this)(k, l) + k - i + 1; + } + + for (int i = 1; i <= l - 1; ++i) { + result *= (*this)(i, l - 1) - (*this)(k, l) + k - i; + } + + for (int i = 1; i <= l; ++i) { + if (i == k) continue; + result /= (*this)(i, l) - (*this)(k, l) + k - i + 1; + result /= (*this)(i, l) - (*this)(k, l) + k - i; + } + + return std::sqrt(-result); +} + +double clebsch::pattern::raising_coeff(int k, int l) const { + double result = 1.0; + + for (int i = 1; i <= l + 1; ++i) { + result *= (*this)(i, l + 1) - (*this)(k, l) + k - i; + } + + for (int i = 1; i <= l - 1; ++i) { + result *= (*this)(i, l - 1) - (*this)(k, l) + k - i - 1; + } + + for (int i = 1; i <= l; ++i) { + if (i == k) continue; + result /= (*this)(i, l) - (*this)(k, l) + k - i; + result /= (*this)(i, l) - (*this)(k, l) + k - i - 1; + } + + return std::sqrt(-result); +} + +// implementation of "decomposition" starts here + +clebsch::decomposition::decomposition(const weight &factor1, const weight &factor2) : + N(factor1.N), factor1(factor1), factor2(factor2) { + assert(factor1.N == factor2.N); + std::vector result; + pattern low(factor1), high(factor1); + weight trial(factor2); + int k = 1, l = N; + + do { + while (k <= N) { + --l; + if (k <= l) { + low(k, l) = std::max(high(k + N - l, N), high(k, l + 1) + trial(l + 1) - trial(l)); + high(k, l) = high(k, l + 1); + if (k > 1 && high(k, l) > high(k - 1, l - 1)) { + high(k, l) = high(k - 1, l - 1); + } + if (l > 1 && k == l && high(k, l) > trial(l - 1) - trial(l)) { + high(k, l) = trial(l - 1) - trial(l); + } + if (low(k, l) > high(k, l)) { + break; + } + trial(l + 1) += high(k, l + 1) - high(k, l); + } else { + trial(l + 1) += high(k, l + 1); + ++k; + l = N; + } + } + + if (k > N) { + result.push_back(trial); + for (int i = 1; i <= N; ++i) { + result.back()(i) -= result.back()(N); + } + } else { + ++l; + } + + while (k != 1 || l != N) { + if (l == N) { + l = --k - 1; + trial(l + 1) -= high(k, l + 1); + } else if (low(k, l) < high(k, l)) { + --high(k, l); + ++trial(l + 1); + break; + } else { + trial(l + 1) -= high(k, l + 1) - high(k, l); + } + ++l; + } + } while (k != 1 || l != N); + + sort(result.begin(), result.end()); + for (std::vector::iterator it = result.begin(); it != result.end(); ++it) { + if (it != result.begin() && *it == weights.back()) { + ++multiplicities.back(); + } else { + weights.push_back(*it); + multiplicities.push_back(1); + } + } +} + +int clebsch::decomposition::size() const { + return weights.size(); +} + +const clebsch::weight &clebsch::decomposition::operator()(int j) const { + return weights[j]; +} + +int clebsch::decomposition::multiplicity(const weight &irrep) const { + assert(irrep.N == N); + std::vector::const_iterator it + = std::lower_bound(weights.begin(), weights.end(), irrep); + + return it != weights.end() && *it == irrep ? multiplicities[it - weights.begin()] : 0; +} + +// implementation of "index_adapter" starts here + +clebsch::index_adapter::index_adapter(const clebsch::decomposition &decomp) : + N(decomp.N), + factor1(decomp.factor1.index()), + factor2(decomp.factor2.index()) { + for (int i = 0, s = decomp.size(); i < s; ++i) { + indices.push_back(decomp(i).index()); + multiplicities.push_back(decomp.multiplicity(decomp(i))); + } +} + +int clebsch::index_adapter::size() const { + return indices.size(); +} + +int clebsch::index_adapter::operator()(int j) const { + return indices[j]; +} + +int clebsch::index_adapter::multiplicity(int irrep) const { + std::vector::const_iterator it = std::lower_bound(indices.begin(), indices.end(), irrep); + + return it != indices.end() && *it == irrep ? multiplicities[it - indices.begin()] : 0; +} + +// implementation of "clebsch" starts here + +void clebsch::coefficients::set(int factor1_state, + int factor2_state, + int multiplicity_index, + int irrep_state, + double value) { + assert(0 <= factor1_state && factor1_state < factor1_dimension); + assert(0 <= factor2_state && factor2_state < factor2_dimension); + assert(0 <= multiplicity_index && multiplicity_index < multiplicity); + assert(0 <= irrep_state && irrep_state < irrep_dimension); + + int coefficient_label[] = { factor1_state, + factor2_state, + multiplicity_index, + irrep_state }; + clzx[std::vector(coefficient_label, coefficient_label + + sizeof coefficient_label / sizeof coefficient_label[0])] = value; +} + +void clebsch::coefficients::highest_weight_normal_form() { + int hws = irrep_dimension - 1; + + // bring CGCs into reduced row echelon form + for (int h = 0, i = 0; h < multiplicity - 1 && i < factor1_dimension; ++i) { + for (int j = 0; h < multiplicity - 1 && j < factor2_dimension; ++j) { + int k0 = h; + + for (int k = h + 1; k < multiplicity; ++k) { + if (fabs((*this)(i, j, k, hws)) > fabs((*this)(i, j, k0, hws))) { + k0 = k; + } + } + + if ((*this)(i, j, k0, hws) < -EPS) { + for (int i2 = i; i2 < factor1_dimension; ++i2) { + for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) { + set(i2, j2, k0, hws, -(*this)(i2, j2, k0, hws)); + } + } + } else if ((*this)(i, j, k0, hws) < EPS) { + continue; + } + + if (k0 != h) { + for (int i2 = i; i2 < factor1_dimension; ++i2) { + for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) { + double x = (*this)(i2, j2, k0, hws); + set(i2, j2, k0, hws, (*this)(i2, j2, h, hws)); + set(i2, j2, h, hws, x); + } + } + } + + for (int k = h + 1; k < multiplicity; ++k) { + for (int i2 = i; i2 < factor1_dimension; ++i2) { + for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) { + set(i2, j2, k, hws, (*this)(i2, j2, k, hws) - (*this)(i2, j2, h, hws) + * (*this)(i, j, k, hws) / (*this)(i, j, h, hws)); + } + } + } + + // next 3 lines not strictly necessary, might improve numerical stability + for (int k = h + 1; k < multiplicity; ++k) { + set(i, j, k, hws, 0.0); + } + + ++h; + } + } + + // Gram-Schmidt orthonormalization + for (int h = 0; h < multiplicity; ++h) { + for (int k = 0; k < h; ++k) { + double overlap = 0.0; + for (int i = 0; i < factor1_dimension; ++i) { + for (int j = 0; j < factor2_dimension; ++j) { + overlap += (*this)(i, j, h, hws) * (*this)(i, j, k, hws); + } + } + + for (int i = 0; i < factor1_dimension; ++i) { + for (int j = 0; j < factor2_dimension; ++j) { + set(i, j, h, hws, (*this)(i, j, h, hws) - overlap * (*this)(i, j, k, hws)); + } + } + } + + double norm = 0.0; + for (int i = 0; i < factor1_dimension; ++i) { + for (int j = 0; j < factor2_dimension; ++j) { + norm += (*this)(i, j, h, hws) * (*this)(i, j, h, hws); + } + } + norm = std::sqrt(norm); + + for (int i = 0; i < factor1_dimension; ++i) { + for (int j = 0; j < factor2_dimension; ++j) { + set(i, j, h, hws, (*this)(i, j, h, hws) / norm); + } + } + } +} + +void clebsch::coefficients::compute_highest_weight_coeffs() { + if (multiplicity == 0) { + return; + } + + std::vector > map_coeff(factor1_dimension, + std::vector(factor2_dimension, -1)); + std::vector > map_states(factor1_dimension, + std::vector(factor2_dimension, -1)); + int n_coeff = 0, n_states = 0; + pattern p(factor1, 0); + + for (int i = 0; i < factor1_dimension; ++i, ++p) { + weight pw(p.get_weight()); + pattern q(factor2, 0); + for (int j = 0; j < factor2_dimension; ++j, ++q) { + if (pw + q.get_weight() == irrep) { + map_coeff[i][j] = n_coeff++; + } + } + } + + if (n_coeff == 1) { + for (int i = 0; i < factor1_dimension; ++i) { + for (int j = 0; j < factor2_dimension; ++j) { + if (map_coeff[i][j] >= 0) { + set(i, j, 0, irrep_dimension - 1, 1.0); + return; + } + } + } + } + + double *hw_system = new double[n_coeff * (factor1_dimension * factor2_dimension)]; + assert(hw_system != NULL); + memset(hw_system, 0, n_coeff * (factor1_dimension * factor2_dimension) * sizeof (double)); + + pattern r(factor1, 0); + for (int i = 0; i < factor1_dimension; ++i, ++r) { + pattern q(factor2, 0); + + for (int j = 0; j < factor2_dimension; ++j, ++q) { + if (map_coeff[i][j] >= 0) { + for (int l = 1; l <= N - 1; ++l) { + for (int k = 1; k <= l; ++k) { + if ((k == 1 || r(k, l) + 1 <= r(k - 1, l - 1)) && r(k, l) + 1 <= r(k, l + 1)) { + ++r(k, l); + int h = r.index(); + --r(k, l); + + if (map_states[h][j] < 0) { + map_states[h][j] = n_states++; + } + + hw_system[n_coeff * map_states[h][j] + map_coeff[i][j]] + += r.raising_coeff(k, l); + } + + if ((k == 1 || q(k, l) + 1 <= q(k - 1, l - 1)) && q(k, l) + 1 <= q(k, l + 1)) { + ++q(k, l); + int h = q.index(); + --q(k, l); + + if (map_states[i][h] < 0) { + map_states[i][h] = n_states++; + } + + + hw_system[n_coeff * map_states[i][h] + map_coeff[i][j]] + += q.raising_coeff(k, l); + } + } + } + } + } + } + + int lwork = -1, info; + double worksize; + + double *singval = new double[std::min(n_coeff, n_states)]; + assert(singval != NULL); + double *singvec = new double[n_coeff * n_coeff]; + assert(singvec != NULL); + + dgesvd_("A", + "N", + &n_coeff, + &n_states, + hw_system, + &n_coeff, + singval, + singvec, + &n_coeff, + NULL, + &n_states, + &worksize, + &lwork, + &info); + assert(info == 0); + + lwork = worksize; + double *work = new double[lwork]; + assert(work != NULL); + + dgesvd_("A", + "N", + &n_coeff, + &n_states, + hw_system, + &n_coeff, + singval, + singvec, + &n_coeff, + NULL, + &n_states, + work, + &lwork, + &info); + assert(info == 0); + + for (int i = 0; i < multiplicity; ++i) { + for (int j = 0; j < factor1_dimension; ++j) { + for (int k = 0; k < factor2_dimension; ++k) { + if (map_coeff[j][k] >= 0) { + double x = singvec[n_coeff * (n_coeff - 1 - i) + map_coeff[j][k]]; + + if (fabs(x) > EPS) { + set(j, k, i, irrep_dimension - 1, x); + } + } + } + } + } + + // uncomment next line to bring highest-weight coefficients into "normal form" + // highest_weight_normal_form(); + + delete[] work; + delete[] singvec; + delete[] singval; + delete[] hw_system; +} + +void clebsch::coefficients::compute_lower_weight_coeffs(int multip_index, + int state, + std::vector &done) { + weight statew(pattern(irrep, state).get_weight()); + pattern p(irrep, 0); + std::vector map_parent(irrep_dimension, -1), + map_multi(irrep_dimension, -1), + which_l(irrep_dimension, -1); + int n_parent = 0, n_multi = 0; + + for (int i = 0; i < irrep_dimension; ++i, ++p) { + weight v(p.get_weight()); + + if (v == statew) { + map_multi[i] = n_multi++; + } else for (int l = 1; l < N; ++l) { + --v(l); + ++v(l + 1); + if (v == statew) { + map_parent[i] = n_parent++; + which_l[i] = l; + if (!done[i]) { + compute_lower_weight_coeffs(multip_index, i, done); + } + break; + } + --v(l + 1); + ++v(l); + } + } + + double *irrep_coeffs = new double[n_parent * n_multi]; + assert(irrep_coeffs != NULL); + memset(irrep_coeffs, 0, n_parent * n_multi * sizeof (double)); + + double *prod_coeffs = new double[n_parent * factor1_dimension * factor2_dimension]; + assert(prod_coeffs != NULL); + memset(prod_coeffs, 0, n_parent * factor1_dimension * factor2_dimension * sizeof (double)); + + std::vector > map_prodstat(factor1_dimension, + std::vector(factor2_dimension, -1)); + int n_prodstat = 0; + + pattern r(irrep, 0); + for (int i = 0; i < irrep_dimension; ++i, ++r) { + if (map_parent[i] >= 0) { + for (int k = 1, l = which_l[i]; k <= l; ++k) { + if (r(k, l) > r(k + 1, l + 1) && (k == l || r(k, l) > r(k, l - 1))) { + --r(k, l); + int h = r.index(); + ++r(k, l); + + irrep_coeffs[n_parent * map_multi[h] + map_parent[i]] += r.lowering_coeff(k, l); + } + } + + pattern q1(factor1, 0); + for (int j1 = 0; j1 < factor1_dimension; ++j1, ++q1) { + pattern q2(factor2, 0); + + for (int j2 = 0; j2 < factor2_dimension; ++j2, ++q2) { + if (std::fabs((*this)(j1, j2, multip_index, i)) > EPS) { + for (int k = 1, l = which_l[i]; k <= l; ++k) { + if (q1(k, l) > q1(k + 1, l + 1) && (k == l || q1(k, l) > q1(k, l - 1))) { + --q1(k, l); + int h = q1.index(); + ++q1(k, l); + + if (map_prodstat[h][j2] < 0) { + map_prodstat[h][j2] = n_prodstat++; + } + + prod_coeffs[n_parent * map_prodstat[h][j2] + map_parent[i]] += + (*this)(j1, j2, multip_index, i) * q1.lowering_coeff(k, l); + } + + if (q2(k, l) > q2(k + 1, l + 1) && (k == l || q2(k, l) > q2(k, l - 1))) { + --q2(k, l); + int h = q2.index(); + ++q2(k, l); + + if (map_prodstat[j1][h] < 0) { + map_prodstat[j1][h] = n_prodstat++; + } + + prod_coeffs[n_parent * map_prodstat[j1][h] + map_parent[i]] += + (*this)(j1, j2, multip_index, i) * q2.lowering_coeff(k, l); + } + } + } + } + } + } + } + + double worksize; + int lwork = -1, info; + + dgels_("N", + &n_parent, + &n_multi, + &n_prodstat, + irrep_coeffs, + &n_parent, + prod_coeffs, + &n_parent, + &worksize, + &lwork, + &info); + assert(info == 0); + + lwork = worksize; + double *work = new double[lwork]; + assert(work != NULL); + + dgels_("N", + &n_parent, + &n_multi, + &n_prodstat, + irrep_coeffs, + &n_parent, + prod_coeffs, + &n_parent, + work, + &lwork, + &info); + assert(info == 0); + + for (int i = 0; i < irrep_dimension; ++i) { + if (map_multi[i] >= 0) { + for (int j = 0; j < factor1_dimension; ++j) { + for (int k = 0; k < factor2_dimension; ++k) { + if (map_prodstat[j][k] >= 0) { + double x = prod_coeffs[n_parent * map_prodstat[j][k] + map_multi[i]]; + + if (fabs(x) > EPS) { + set(j, k, multip_index, i, x); + } + } + } + } + + done[i] = true; + } + } + + delete[] work; + delete[] prod_coeffs; + delete[] irrep_coeffs; +} + +clebsch::coefficients::coefficients(const weight &irrep, const weight &factor1, const weight &factor2) : + N(irrep.N), + factor1(factor1), + factor2(factor2), + irrep(irrep), + factor1_dimension(factor1.dimension()), + factor2_dimension(factor2.dimension()), + irrep_dimension(irrep.dimension()), + multiplicity(decomposition(factor1, factor2).multiplicity(irrep)) { + assert(factor1.N == irrep.N); + assert(factor2.N == irrep.N); + + compute_highest_weight_coeffs(); + + for (int i = 0; i < multiplicity; ++i) { + std::vector done(irrep_dimension, 0); + done[irrep_dimension - 1] = true; + for (int j = irrep_dimension - 1; j >= 0; --j) { + if (!done[j]) { + compute_lower_weight_coeffs(i, j, done); + } + } + } +} + +double clebsch::coefficients::operator()(int factor1_state, + int factor2_state, + int multiplicity_index, + int irrep_state) const { + + assert(0 <= factor1_state && factor1_state < factor1_dimension); + assert(0 <= factor2_state && factor2_state < factor2_dimension); + assert(0 <= multiplicity_index && multiplicity_index < multiplicity); + assert(0 <= irrep_state && irrep_state < irrep_dimension); + + int coefficient_label[] = { factor1_state, + factor2_state, + multiplicity_index, + irrep_state }; + std::map, double>::const_iterator it( + clzx.find(std::vector(coefficient_label, coefficient_label + + sizeof coefficient_label / sizeof coefficient_label[0]))); + + return it != clzx.end() ? it->second : 0.0; +} + diff --git a/clebsch.h b/clebsch.h new file mode 100644 index 0000000..21fe4ce --- /dev/null +++ b/clebsch.h @@ -0,0 +1,236 @@ +#ifndef CLEBSCH_H +#define CLEBSCH_H + +//#define NDEBUG + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Declaration of LAPACK subroutines +// Make sure the data types match your version of LAPACK + +extern "C" void dgesvd_(char const* JOBU, + char const* JOBVT, + int const* M, + int const* N, + double* A, + int const* LDA, + double* S, + double* U, + int const* LDU, + double* VT, + int const* LDVT, + double* WORK, + int const* LWORK, + int *INFO); + +extern "C" void dgels_(char const* TRANS, + int const* M, + int const* N, + int const* NRHS, + double* A, + int const* LDA, + double* B, + int const* LDB, + double* WORK, + int const* LWORK, + int *INFO); + +namespace clebsch { + const double EPS = 1e-12; + + // binomial coefficients + class binomial_t { + std::vector cache; + int N; + + public: + int operator()(int n, int k); + }; + + // Eq. (19) and (25) + class weight { + std::vector elem; + + public: + // the N in "SU(N)" + const int N; + + // create a non-initialized weight + weight(int N); + + // create irrep weight of given index + // Eq. (C2) + weight(int N, int index); + + // assign from another instance + clebsch::weight &operator=(const clebsch::weight &w); + + // access elements of this weight (k = 1, ..., N) + int &operator()(int k); + const int &operator()(int k) const; + + // compare weights + // Eq. (C1) + bool operator<(const weight &w) const; + bool operator==(const weight &w) const; + + // element-wise sum of weights + clebsch::weight operator+(const weight &w) const; + + // returns the index of this irrep weight (index = 0, 1, ...) + // Eq. (C2) + int index() const; + + // returns the dimension of this irrep weight + // Eq. (22) + long long dimension() const; + }; + + // Eq. (20) + class pattern { + std::vector elem; + + public: + // the N in "SU(N)" + const int N; + + // copy constructor + pattern(const pattern &pat); + + // create pattern of given index from irrep weight + // Eq. (C7) + pattern(const weight &irrep, int index = 0); + + // access elements of this pattern (l = 1, ..., N; k = 1, ..., l) + int &operator()(int k, int l); + const int &operator()(int k, int l) const; + + // find succeeding/preceding pattern, return false if not possible + // Eq. (C9) + bool operator++(); + bool operator--(); + + // returns the pattern index (index = 0, ..., dimension - 1) + // Eq. (C7) + int index() const; + + // returns the pattern weight + // Eq. (25) + clebsch::weight get_weight() const; + + // returns matrix element of lowering operator J^(l)_- + // between this pattern minus M^(k,l) and this pattern + // (l = 1, ..., N; k = 1, ..., l) + // Eq. (28) + double lowering_coeff(int k, int l) const; + + // returns matrix element of raising operator J^(l)_+ + // between this pattern plus M^(k,l) and this pattern + // (l = 1, ..., N; k = 1, ..., l) + // Eq. (29) + double raising_coeff(int k, int l) const; + }; + + class decomposition { + std::vector weights; + std::vector multiplicities; + + public: + // the N in "SU(N)" + const int N; + + // save given irreps for later use + const weight factor1, factor2; + + // construct the decomposition of factor1 times factor2 into irreps + // Eq. (31) + decomposition(const weight &factor1, const weight &factor2); + + // return the number of occurring irreps + int size() const; + + // access the occurring irreps + // j = 0, ..., size() - 1 + const clebsch::weight &operator()(int j) const; + + // return the outer multiplicity of irrep in this decomposition + int multiplicity(const weight &irrep) const; + }; + + class index_adapter { + std::vector indices; + std::vector multiplicities; + + public: + // the N in "SU(N)" + const int N; + + // save given irreps for later use + const int factor1, factor2; + + // construct this index_adapter from a given decomposition + index_adapter(const clebsch::decomposition &decomp); + + // return the number of occurring irreps + int size() const; + + // access the occurring irreps + int operator()(int j) const; + + // return the outer multiplicity of irrep in this decomposition + int multiplicity(int irrep) const; + }; + + class coefficients { + std::map, double> clzx; + + // access Clebsch-Gordan coefficients in convenient manner + void set(int factor1_state, + int factor2_state, + int multiplicity_index, + int irrep_state, + double value); + + // internal functions, doing most of the work + void highest_weight_normal_form(); // Eq. (37) + void compute_highest_weight_coeffs(); // Eq. (36) + void compute_lower_weight_coeffs(int multip_index, int state, std::vector &done); // Eq. (40) + + public: + // the N in "SU(N)" + const int N; + + // save irreps and their dimensions for later use + const weight factor1, factor2, irrep; + const int factor1_dimension, factor2_dimension, irrep_dimension; + + // outer multiplicity of irrep in this decomposition + const int multiplicity; + + // construct all Clebsch-Gordan coefficients of this decomposition + coefficients(const weight &irrep, const weight &factor1, const weight &factor2); + + // access Clebsch-Gordan coefficients (read-only) + // multiplicity_index = 0, ..., multiplicity - 1 + // factor1_state = 0, ..., factor1_dimension - 1 + // factor2_state = 0, ..., factor2_dimension - 1 + // irrep_state = 0, ..., irrep_dimension + double operator()(int factor1_state, + int factor2_state, + int multiplicity_index, + int irrep_state) const; + }; +}; + +#endif diff --git a/clebschgordan.cpp b/clebschgordan.cpp new file mode 100644 index 0000000..81a37bd --- /dev/null +++ b/clebschgordan.cpp @@ -0,0 +1,251 @@ +#include + +#include "clebsch.h" + +using namespace std; + +int main() { + while (true) { + int choice, N; + + cout << "What would you like to do?" << endl; + cout << "1) Translate an i-weight S to its index P(S)" << endl; + cout << "2) Recover an i-weight S from its index P(S)" << endl; + cout << "3) Translate a pattern M to its index Q(M)" << endl; + cout << "4) Recover a pattern M from its index Q(M)" << endl; + cout << "5) Calculate Clebsch-Gordan coefficients for S x S' -> S''" << endl; + cout << "6) Calculate all Glebsch-Gordan coefficients for S x S'" << endl; + cout << "0) Quit" << endl; + + do { + cin >> choice; + } while (choice < 0 || choice > 6); + + if (choice == 0) { + break; + } + + cout << "N (e.g. 3): "; + cin >> N; + + switch (choice) { + case 1: { + clebsch::weight S(N); + cout << "Irrep S: "; + for (int k = 1; k <= N; ++k) { + cin >> S(k); + } + cout << S.index() << endl; + break; + } + case 2: { + int P; + cout << "Index: "; + cin >> P; + clebsch::weight S(N, P); + cout << "I-weight:"; + for (int k = 1; k <= N; ++k) { + cout << ' ' << S(k); + } + cout << endl; + break; + } + case 3: { + clebsch::pattern M(N); + for (int l = N; l >= 1; --l) { + cout << "Row l = " << l << ": "; + for (int k = 1; k <= l; ++k) { + cin >> M(k, l); + } + } + cout << "Index: " << M.index() + 1 << endl; + break; + } + case 4: { + clebsch::weight S(N); + cout << "Irrep S: "; + for (int i = 1; i <= N; ++i) { + cin >> S(i); + } + + int Q; + cout << "Index (1..dim(S)): "; + cin >> Q; + clebsch::pattern M(S, Q - 1); + for (int l = N; l >= 1; --l) { + for (int k = 1; k <= l; ++k) { + cout << M(k, l) << '\t'; + } + cout << endl; + } + break; + } + case 5: { + clebsch::weight S(N); + cout << "Irrep S (e.g."; + for (int k = N - 1; k >= 0; --k) { + cout << ' ' << k; + } + cout << "): "; + for (int k = 1; k <= N; ++k) { + cin >> S(k); + } + + clebsch::weight Sprime(N); + cout << "Irrep S' (e.g."; + for (int k = N - 1; k >= 0; --k) { + cout << ' ' << k; + } + cout << "): "; + for (int k = 1; k <= N; ++k) { + cin >> Sprime(k); + } + + clebsch::decomposition decomp(S, Sprime); + cout << "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl; + cout << "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl; + for (int i = 0; i < decomp.size(); ++i) { + cout << "[" << decomp(i).index() << "] "; + for (int k = 1; k <= N; ++k) { + cout << decomp(i)(k) << ' '; + } + cout << '(' << decomp.multiplicity(decomp(i)) << ") {" + << decomp(i).dimension() << "}" << endl;; + } + + clebsch::weight Sdoubleprime(N); + for (bool b = true; b; ) { + cout << "Irrep S'': "; + for (int k = 1; k <= N; ++k) { + cin >> Sdoubleprime(k); + } + for (int i = 0; i < decomp.size(); ++i) { + if (decomp(i) == Sdoubleprime) { + b = false; + break; + } + } + if (b) { + cout << "S'' does not occur in the decomposition" << endl; + } + } + + int alpha; + while (true) { + cout << "Outer multiplicity index: "; + cin >> alpha; + if (1 <= alpha && alpha <= decomp.multiplicity(Sdoubleprime)) { + break; + } + cout << "S'' does not have this multiplicity" << endl; + } + + string file_name; + cout << "Enter file name to write to file (leave blank for screen output): "; + cin.ignore(1234, '\n'); + getline(cin, file_name); + + const clebsch::coefficients C(Sdoubleprime, S, Sprime); + int dimS = S.dimension(), + dimSprime = Sprime.dimension(), + dimSdoubleprime = Sdoubleprime.dimension(); + + ofstream os(file_name.c_str()); + (file_name.empty() ? cout : os).setf(ios::fixed); + (file_name.empty() ? cout : os).precision(15); + (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'', alpha" << endl; + (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl; + for (int i = 0; i < dimSdoubleprime; ++i) { + for (int j = 0; j < dimS; ++j) { + for (int k = 0; k < dimSprime; ++k) { + double x = double(C(j, k, alpha - 1, i)); + + if (fabs(x) > clebsch::EPS) { + (file_name.empty() ? cout : os) << j + 1 << '\t' + << k + 1 << '\t' << i + 1 << '\t' << x << endl; + } + } + } + } + + break; + } + case 6: { + clebsch::weight S(N); + cout << "Irrep S (e.g."; + for (int k = N - 1; k >= 0; --k) { + cout << ' ' << k; + } + cout << "): "; + for (int k = 1; k <= N; ++k) { + cin >> S(k); + } + + clebsch::weight Sprime(N); + cout << "Irrep S' (e.g."; + for (int k = N - 1; k >= 0; --k) { + cout << ' ' << k; + } + cout << "): "; + for (int k = 1; k <= N; ++k) { + cin >> Sprime(k); + } + + string file_name; + cout << "Enter file name to write to file (leave blank for screen output): "; + cin.ignore(1234, '\n'); + getline(cin, file_name); + + ofstream os(file_name.c_str()); + (file_name.empty() ? cout : os).setf(ios::fixed); + (file_name.empty() ? cout : os).precision(15); + + clebsch::decomposition decomp(S, Sprime); + (file_name.empty() ? cout : os) << + "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl; + (file_name.empty() ? cout : os) << + "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl; + for (int i = 0; i < decomp.size(); ++i) { + (file_name.empty() ? cout : os) << "[" << decomp(i).index() << "] "; + for (int k = 1; k <= N; ++k) { + (file_name.empty() ? cout : os) << decomp(i)(k) << ' '; + } + (file_name.empty() ? cout : os) << '(' << decomp.multiplicity(decomp(i)) << ") {" + << decomp(i).dimension() << "}" << endl;; + } + + for (int i = 0; i < decomp.size(); ++i) { + const clebsch::coefficients C(decomp(i),S, Sprime); + int dimS = S.dimension(), + dimSprime = Sprime.dimension(), + dimSdoubleprime = decomp(i).dimension(); + + for (int m = 0; m < C.multiplicity; ++m) { + (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'' = ("; + for (int j = 1; j <= N; ++j) cout << decomp(i)(j) << (j < N ? ' ' : ')'); + (file_name.empty() ? cout : os) << ", alpha = " << m + 1 << endl; + (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl; + for (int i = 0; i < dimSdoubleprime; ++i) { + for (int j = 0; j < dimS; ++j) { + for (int k = 0; k < dimSprime; ++k) { + double x = double(C(j, k, m, i)); + + if (fabs(x) > clebsch::EPS) { + (file_name.empty() ? cout : os) << j + 1<< '\t' + << k + 1 << '\t' << i + 1 << '\t' << x << endl; + } + } + } + } + + (file_name.empty() ? cout : os) << endl; + } + } + + break; + } + } + } + + return 0; +} diff --git a/su2clebsch.cpp b/su2clebsch.cpp new file mode 100644 index 0000000..49a8047 --- /dev/null +++ b/su2clebsch.cpp @@ -0,0 +1,60 @@ +#include + +#include "clebsch.h" + +#include "su2clebsch.hpp" + +using namespace std; + +namespace su2clebsch { + + double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, + int irrep_2j, int irrep_2m) + { + 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; + for (int i = 0; i < decomp.size(); ++i) { + if (decomp(i) == irrep_S) { + b = false; + break; + } + } + if (b) { + cerr << "WARNING: 2j=" << irrep_2j + << " does not occour in decomposition!" << endl; + return 0.0; + } + + const clebsch::coefficients C(irrep_S, factor1_S, factor2_S); + + return C( factor1_n, factor2_n, 0, irrep_n ); + } + +}; + diff --git a/su2clebsch.hpp b/su2clebsch.hpp new file mode 100644 index 0000000..689c0a1 --- /dev/null +++ b/su2clebsch.hpp @@ -0,0 +1,11 @@ +#ifndef SU2CLEBSCH_HPP +#define SU2CLEBSCH_HPP + +namespace su2clebsch { + + double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m, + int irrep_2j, int irrep_2m); + +}; + +#endif