294 lines
11 KiB
C++
294 lines
11 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at the mozilla.org home page
|
|
|
|
#ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
|
|
#define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
|
|
|
|
namespace Eigen {
|
|
|
|
class DynamicSGroup
|
|
{
|
|
public:
|
|
inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) { m_elements.push_back(ge(Generator(0, 0, 0))); }
|
|
inline DynamicSGroup(const DynamicSGroup& o) : m_numIndices(o.m_numIndices), m_elements(o.m_elements), m_generators(o.m_generators), m_globalFlags(o.m_globalFlags) { }
|
|
inline DynamicSGroup(DynamicSGroup&& o) : m_numIndices(o.m_numIndices), m_elements(), m_generators(o.m_generators), m_globalFlags(o.m_globalFlags) { std::swap(m_elements, o.m_elements); }
|
|
inline DynamicSGroup& operator=(const DynamicSGroup& o) { m_numIndices = o.m_numIndices; m_elements = o.m_elements; m_generators = o.m_generators; m_globalFlags = o.m_globalFlags; return *this; }
|
|
inline DynamicSGroup& operator=(DynamicSGroup&& o) { m_numIndices = o.m_numIndices; std::swap(m_elements, o.m_elements); m_generators = o.m_generators; m_globalFlags = o.m_globalFlags; return *this; }
|
|
|
|
void add(int one, int two, int flags = 0);
|
|
|
|
template<typename Gen_>
|
|
inline void add(Gen_) { add(Gen_::One, Gen_::Two, Gen_::Flags); }
|
|
inline void addSymmetry(int one, int two) { add(one, two, 0); }
|
|
inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
|
|
inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
|
|
inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); }
|
|
|
|
template<typename Op, typename RV, typename Index, std::size_t N, typename... Args>
|
|
inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const
|
|
{
|
|
eigen_assert(N >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
|
|
for (std::size_t i = 0; i < size(); i++)
|
|
initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags, initial, std::forward<Args>(args)...);
|
|
return initial;
|
|
}
|
|
|
|
template<typename Op, typename RV, typename Index, typename... Args>
|
|
inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const
|
|
{
|
|
eigen_assert(idx.size() >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
|
|
for (std::size_t i = 0; i < size(); i++)
|
|
initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
|
|
return initial;
|
|
}
|
|
|
|
inline int globalFlags() const { return m_globalFlags; }
|
|
inline std::size_t size() const { return m_elements.size(); }
|
|
|
|
template<typename Tensor_, typename... IndexTypes>
|
|
inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, typename Tensor_::Index firstIndex, IndexTypes... otherIndices) const
|
|
{
|
|
static_assert(sizeof...(otherIndices) + 1 == Tensor_::NumIndices, "Number of indices used to access a tensor coefficient must be equal to the rank of the tensor.");
|
|
return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
|
|
}
|
|
|
|
template<typename Tensor_>
|
|
inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const
|
|
{
|
|
return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
|
|
}
|
|
private:
|
|
struct GroupElement {
|
|
std::vector<int> representation;
|
|
int flags;
|
|
bool isId() const
|
|
{
|
|
for (std::size_t i = 0; i < representation.size(); i++)
|
|
if (i != (size_t)representation[i])
|
|
return false;
|
|
return true;
|
|
}
|
|
};
|
|
struct Generator {
|
|
int one;
|
|
int two;
|
|
int flags;
|
|
constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
|
|
};
|
|
|
|
std::size_t m_numIndices;
|
|
std::vector<GroupElement> m_elements;
|
|
std::vector<Generator> m_generators;
|
|
int m_globalFlags;
|
|
|
|
template<typename Index, std::size_t N, int... n>
|
|
inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx, internal::numeric_list<int, n...>) const
|
|
{
|
|
return std::array<Index, N>{{ idx[n >= m_numIndices ? n : m_elements[which].representation[n]]... }};
|
|
}
|
|
|
|
template<typename Index>
|
|
inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const
|
|
{
|
|
std::vector<Index> result;
|
|
result.reserve(idx.size());
|
|
for (auto k : m_elements[which].representation)
|
|
result.push_back(idx[k]);
|
|
for (std::size_t i = m_numIndices; i < idx.size(); i++)
|
|
result.push_back(idx[i]);
|
|
return result;
|
|
}
|
|
|
|
inline GroupElement ge(Generator const& g) const
|
|
{
|
|
GroupElement result;
|
|
result.representation.reserve(m_numIndices);
|
|
result.flags = g.flags;
|
|
for (std::size_t k = 0; k < m_numIndices; k++) {
|
|
if (k == (std::size_t)g.one)
|
|
result.representation.push_back(g.two);
|
|
else if (k == (std::size_t)g.two)
|
|
result.representation.push_back(g.one);
|
|
else
|
|
result.representation.push_back(int(k));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
GroupElement mul(GroupElement, GroupElement) const;
|
|
inline GroupElement mul(Generator g1, GroupElement g2) const
|
|
{
|
|
return mul(ge(g1), g2);
|
|
}
|
|
|
|
inline GroupElement mul(GroupElement g1, Generator g2) const
|
|
{
|
|
return mul(g1, ge(g2));
|
|
}
|
|
|
|
inline GroupElement mul(Generator g1, Generator g2) const
|
|
{
|
|
return mul(ge(g1), ge(g2));
|
|
}
|
|
|
|
inline int findElement(GroupElement e) const
|
|
{
|
|
for (auto ee : m_elements) {
|
|
if (ee.representation == e.representation)
|
|
return ee.flags ^ e.flags;
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
void updateGlobalFlags(int flagDiffOfSameGenerator);
|
|
};
|
|
|
|
// dynamic symmetry group that auto-adds the template parameters in the constructor
|
|
template<typename... Gen>
|
|
class DynamicSGroupFromTemplateArgs : public DynamicSGroup
|
|
{
|
|
public:
|
|
inline DynamicSGroupFromTemplateArgs() : DynamicSGroup()
|
|
{
|
|
add_all(internal::type_list<Gen...>());
|
|
}
|
|
inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) { }
|
|
inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) { }
|
|
inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) { DynamicSGroup::operator=(o); return *this; }
|
|
inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) { DynamicSGroup::operator=(o); return *this; }
|
|
|
|
private:
|
|
template<typename Gen1, typename... GenNext>
|
|
inline void add_all(internal::type_list<Gen1, GenNext...>)
|
|
{
|
|
add(Gen1());
|
|
add_all(internal::type_list<GenNext...>());
|
|
}
|
|
|
|
inline void add_all(internal::type_list<>)
|
|
{
|
|
}
|
|
};
|
|
|
|
inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const
|
|
{
|
|
eigen_internal_assert(g1.representation.size() == m_numIndices);
|
|
eigen_internal_assert(g2.representation.size() == m_numIndices);
|
|
|
|
GroupElement result;
|
|
result.representation.reserve(m_numIndices);
|
|
for (std::size_t i = 0; i < m_numIndices; i++) {
|
|
int v = g2.representation[g1.representation[i]];
|
|
eigen_assert(v >= 0);
|
|
result.representation.push_back(v);
|
|
}
|
|
result.flags = g1.flags ^ g2.flags;
|
|
return result;
|
|
}
|
|
|
|
inline void DynamicSGroup::add(int one, int two, int flags)
|
|
{
|
|
eigen_assert(one >= 0);
|
|
eigen_assert(two >= 0);
|
|
eigen_assert(one != two);
|
|
|
|
if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) {
|
|
std::size_t newNumIndices = (one > two) ? one : two + 1;
|
|
for (auto& gelem : m_elements) {
|
|
gelem.representation.reserve(newNumIndices);
|
|
for (std::size_t i = m_numIndices; i < newNumIndices; i++)
|
|
gelem.representation.push_back(i);
|
|
}
|
|
m_numIndices = newNumIndices;
|
|
}
|
|
|
|
Generator g{one, two, flags};
|
|
GroupElement e = ge(g);
|
|
|
|
/* special case for first generator */
|
|
if (m_elements.size() == 1) {
|
|
while (!e.isId()) {
|
|
m_elements.push_back(e);
|
|
e = mul(e, g);
|
|
}
|
|
|
|
if (e.flags > 0)
|
|
updateGlobalFlags(e.flags);
|
|
|
|
// only add in case we didn't have identity
|
|
if (m_elements.size() > 1)
|
|
m_generators.push_back(g);
|
|
return;
|
|
}
|
|
|
|
int p = findElement(e);
|
|
if (p >= 0) {
|
|
updateGlobalFlags(p);
|
|
return;
|
|
}
|
|
|
|
std::size_t coset_order = m_elements.size();
|
|
m_elements.push_back(e);
|
|
for (std::size_t i = 1; i < coset_order; i++)
|
|
m_elements.push_back(mul(m_elements[i], e));
|
|
m_generators.push_back(g);
|
|
|
|
std::size_t coset_rep = coset_order;
|
|
do {
|
|
for (auto g : m_generators) {
|
|
e = mul(m_elements[coset_rep], g);
|
|
p = findElement(e);
|
|
if (p < 0) {
|
|
// element not yet in group
|
|
m_elements.push_back(e);
|
|
for (std::size_t i = 1; i < coset_order; i++)
|
|
m_elements.push_back(mul(m_elements[i], e));
|
|
} else if (p > 0) {
|
|
updateGlobalFlags(p);
|
|
}
|
|
}
|
|
coset_rep += coset_order;
|
|
} while (coset_rep < m_elements.size());
|
|
}
|
|
|
|
inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator)
|
|
{
|
|
switch (flagDiffOfSameGenerator) {
|
|
case 0:
|
|
default:
|
|
// nothing happened
|
|
break;
|
|
case NegationFlag:
|
|
// every element is it's own negative => whole tensor is zero
|
|
m_globalFlags |= GlobalZeroFlag;
|
|
break;
|
|
case ConjugationFlag:
|
|
// every element is it's own conjugate => whole tensor is real
|
|
m_globalFlags |= GlobalRealFlag;
|
|
break;
|
|
case (NegationFlag | ConjugationFlag):
|
|
// every element is it's own negative conjugate => whole tensor is imaginary
|
|
m_globalFlags |= GlobalImagFlag;
|
|
break;
|
|
/* NOTE:
|
|
* since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator
|
|
* causes the tensor to be real and the next one to be imaginary, this will
|
|
* trivially give the correct result
|
|
*/
|
|
}
|
|
}
|
|
|
|
} // end namespace Eigen
|
|
|
|
#endif // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
|
|
|
|
/*
|
|
* kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
|
|
*/
|