Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/eigen3/unsupported/Eigen/CXX11/src/TensorSymmetry/DynamicSymmetry.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
0011 #define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
0012 
0013 namespace Eigen {
0014 
0015 class DynamicSGroup
0016 {
0017   public:
0018     inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) { m_elements.push_back(ge(Generator(0, 0, 0))); }
0019     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) { }
0020     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); }
0021     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; }
0022     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; }
0023 
0024     void add(int one, int two, int flags = 0);
0025 
0026     template<typename Gen_>
0027     inline void add(Gen_) { add(Gen_::One, Gen_::Two, Gen_::Flags); }
0028     inline void addSymmetry(int one, int two) { add(one, two, 0); }
0029     inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
0030     inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
0031     inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); }
0032 
0033     template<typename Op, typename RV, typename Index, std::size_t N, typename... Args>
0034     inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const
0035     {
0036       eigen_assert(N >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
0037       for (std::size_t i = 0; i < size(); i++)
0038         initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags, initial, std::forward<Args>(args)...);
0039       return initial;
0040     }
0041 
0042     template<typename Op, typename RV, typename Index, typename... Args>
0043     inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const
0044     {
0045       eigen_assert(idx.size() >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
0046       for (std::size_t i = 0; i < size(); i++)
0047         initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
0048       return initial;
0049     }
0050 
0051     inline int globalFlags() const { return m_globalFlags; }
0052     inline std::size_t size() const { return m_elements.size(); }
0053 
0054     template<typename Tensor_, typename... IndexTypes>
0055     inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, typename Tensor_::Index firstIndex, IndexTypes... otherIndices) const
0056     {
0057       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.");
0058       return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
0059     }
0060 
0061     template<typename Tensor_>
0062     inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const
0063     {
0064       return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
0065     }
0066   private:
0067     struct GroupElement {
0068       std::vector<int> representation;
0069       int flags;
0070       bool isId() const
0071       {
0072         for (std::size_t i = 0; i < representation.size(); i++)
0073           if (i != (size_t)representation[i])
0074             return false;
0075         return true;
0076       }
0077     };
0078     struct Generator {
0079       int one;
0080       int two;
0081       int flags;
0082       constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
0083     };
0084 
0085     std::size_t m_numIndices;
0086     std::vector<GroupElement> m_elements;
0087     std::vector<Generator> m_generators;
0088     int m_globalFlags;
0089 
0090     template<typename Index, std::size_t N, int... n>
0091     inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx, internal::numeric_list<int, n...>) const
0092     {
0093       return std::array<Index, N>{{ idx[n >= m_numIndices ? n : m_elements[which].representation[n]]... }};
0094     }
0095 
0096     template<typename Index>
0097     inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const
0098     {
0099       std::vector<Index> result;
0100       result.reserve(idx.size());
0101       for (auto k : m_elements[which].representation)
0102         result.push_back(idx[k]);
0103       for (std::size_t i = m_numIndices; i < idx.size(); i++)
0104         result.push_back(idx[i]);
0105       return result;
0106     }
0107 
0108     inline GroupElement ge(Generator const& g) const
0109     {
0110       GroupElement result;
0111       result.representation.reserve(m_numIndices);
0112       result.flags = g.flags;
0113       for (std::size_t k = 0; k < m_numIndices; k++) {
0114         if (k == (std::size_t)g.one)
0115           result.representation.push_back(g.two);
0116         else if (k == (std::size_t)g.two)
0117           result.representation.push_back(g.one);
0118         else
0119           result.representation.push_back(int(k));
0120       }
0121       return result;
0122     }
0123 
0124     GroupElement mul(GroupElement, GroupElement) const;
0125     inline GroupElement mul(Generator g1, GroupElement g2) const
0126     {
0127       return mul(ge(g1), g2);
0128     }
0129 
0130     inline GroupElement mul(GroupElement g1, Generator g2) const
0131     {
0132       return mul(g1, ge(g2));
0133     }
0134 
0135     inline GroupElement mul(Generator g1, Generator g2) const
0136     {
0137       return mul(ge(g1), ge(g2));
0138     }
0139 
0140     inline int findElement(GroupElement e) const
0141     {
0142       for (auto ee : m_elements) {
0143         if (ee.representation == e.representation)
0144           return ee.flags ^ e.flags;
0145       }
0146       return -1;
0147     }
0148 
0149     void updateGlobalFlags(int flagDiffOfSameGenerator);
0150 };
0151 
0152 // dynamic symmetry group that auto-adds the template parameters in the constructor
0153 template<typename... Gen>
0154 class DynamicSGroupFromTemplateArgs : public DynamicSGroup
0155 {
0156   public:
0157     inline DynamicSGroupFromTemplateArgs() : DynamicSGroup()
0158     {
0159       add_all(internal::type_list<Gen...>());
0160     }
0161     inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) { }
0162     inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) { }
0163     inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) { DynamicSGroup::operator=(o); return *this; }
0164     inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) { DynamicSGroup::operator=(o); return *this; }
0165   
0166   private:
0167     template<typename Gen1, typename... GenNext>
0168     inline void add_all(internal::type_list<Gen1, GenNext...>)
0169     {
0170       add(Gen1());
0171       add_all(internal::type_list<GenNext...>());
0172     }
0173 
0174     inline void add_all(internal::type_list<>)
0175     {
0176     }
0177 };
0178 
0179 inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const
0180 {
0181   eigen_internal_assert(g1.representation.size() == m_numIndices);
0182   eigen_internal_assert(g2.representation.size() == m_numIndices);
0183 
0184   GroupElement result;
0185   result.representation.reserve(m_numIndices);
0186   for (std::size_t i = 0; i < m_numIndices; i++) {
0187     int v = g2.representation[g1.representation[i]];
0188     eigen_assert(v >= 0);
0189     result.representation.push_back(v);
0190   }
0191   result.flags = g1.flags ^ g2.flags;
0192   return result;
0193 }
0194 
0195 inline void DynamicSGroup::add(int one, int two, int flags)
0196 {
0197   eigen_assert(one >= 0);
0198   eigen_assert(two >= 0);
0199   eigen_assert(one != two);
0200 
0201   if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) {
0202     std::size_t newNumIndices = (one > two) ? one : two + 1;
0203     for (auto& gelem : m_elements) {
0204       gelem.representation.reserve(newNumIndices);
0205       for (std::size_t i = m_numIndices; i < newNumIndices; i++)
0206         gelem.representation.push_back(i);
0207     }
0208     m_numIndices = newNumIndices;
0209   }
0210 
0211   Generator g{one, two, flags};
0212   GroupElement e = ge(g);
0213 
0214   /* special case for first generator */
0215   if (m_elements.size() == 1) {
0216     while (!e.isId()) {
0217       m_elements.push_back(e);
0218       e = mul(e, g);
0219     }
0220 
0221     if (e.flags > 0)
0222       updateGlobalFlags(e.flags);
0223 
0224     // only add in case we didn't have identity
0225     if (m_elements.size() > 1)
0226       m_generators.push_back(g);
0227     return;
0228   }
0229 
0230   int p = findElement(e);
0231   if (p >= 0) {
0232     updateGlobalFlags(p);
0233     return;
0234   }
0235 
0236   std::size_t coset_order = m_elements.size();
0237   m_elements.push_back(e);
0238   for (std::size_t i = 1; i < coset_order; i++)
0239     m_elements.push_back(mul(m_elements[i], e));
0240   m_generators.push_back(g);
0241 
0242   std::size_t coset_rep = coset_order;
0243   do {
0244     for (auto g : m_generators) {
0245       e = mul(m_elements[coset_rep], g);
0246       p = findElement(e);
0247       if (p < 0) {
0248         // element not yet in group
0249         m_elements.push_back(e);
0250         for (std::size_t i = 1; i < coset_order; i++)
0251           m_elements.push_back(mul(m_elements[i], e));
0252       } else if (p > 0) {
0253         updateGlobalFlags(p);
0254       }
0255     }
0256     coset_rep += coset_order;
0257   } while (coset_rep < m_elements.size());
0258 }
0259 
0260 inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator)
0261 {
0262     switch (flagDiffOfSameGenerator) {
0263       case 0:
0264       default:
0265         // nothing happened
0266         break;
0267       case NegationFlag:
0268         // every element is it's own negative => whole tensor is zero
0269         m_globalFlags |= GlobalZeroFlag;
0270         break;
0271       case ConjugationFlag:
0272         // every element is it's own conjugate => whole tensor is real
0273         m_globalFlags |= GlobalRealFlag;
0274         break;
0275       case (NegationFlag | ConjugationFlag):
0276         // every element is it's own negative conjugate => whole tensor is imaginary
0277         m_globalFlags |= GlobalImagFlag;
0278         break;
0279       /* NOTE:
0280        *   since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator
0281        *   causes the tensor to be real and the next one to be imaginary, this will
0282        *   trivially give the correct result
0283        */
0284     }
0285 }
0286 
0287 } // end namespace Eigen
0288 
0289 #endif // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
0290 
0291 /*
0292  * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
0293  */