![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2013 Christian Seiler <christian@iwakd.de> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H 00011 #define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H 00012 00013 namespace Eigen { 00014 00015 class DynamicSGroup 00016 { 00017 public: 00018 inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) { m_elements.push_back(ge(Generator(0, 0, 0))); } 00019 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) { } 00020 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); } 00021 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; } 00022 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; } 00023 00024 void add(int one, int two, int flags = 0); 00025 00026 template<typename Gen_> 00027 inline void add(Gen_) { add(Gen_::One, Gen_::Two, Gen_::Flags); } 00028 inline void addSymmetry(int one, int two) { add(one, two, 0); } 00029 inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); } 00030 inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); } 00031 inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); } 00032 00033 template<typename Op, typename RV, typename Index, std::size_t N, typename... Args> 00034 inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const 00035 { 00036 eigen_assert(N >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices."); 00037 for (std::size_t i = 0; i < size(); i++) 00038 initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags, initial, std::forward<Args>(args)...); 00039 return initial; 00040 } 00041 00042 template<typename Op, typename RV, typename Index, typename... Args> 00043 inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const 00044 { 00045 eigen_assert(idx.size() >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices."); 00046 for (std::size_t i = 0; i < size(); i++) 00047 initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...); 00048 return initial; 00049 } 00050 00051 inline int globalFlags() const { return m_globalFlags; } 00052 inline std::size_t size() const { return m_elements.size(); } 00053 00054 template<typename Tensor_, typename... IndexTypes> 00055 inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, typename Tensor_::Index firstIndex, IndexTypes... otherIndices) const 00056 { 00057 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."); 00058 return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}}); 00059 } 00060 00061 template<typename Tensor_> 00062 inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const 00063 { 00064 return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices); 00065 } 00066 private: 00067 struct GroupElement { 00068 std::vector<int> representation; 00069 int flags; 00070 bool isId() const 00071 { 00072 for (std::size_t i = 0; i < representation.size(); i++) 00073 if (i != (size_t)representation[i]) 00074 return false; 00075 return true; 00076 } 00077 }; 00078 struct Generator { 00079 int one; 00080 int two; 00081 int flags; 00082 constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {} 00083 }; 00084 00085 std::size_t m_numIndices; 00086 std::vector<GroupElement> m_elements; 00087 std::vector<Generator> m_generators; 00088 int m_globalFlags; 00089 00090 template<typename Index, std::size_t N, int... n> 00091 inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx, internal::numeric_list<int, n...>) const 00092 { 00093 return std::array<Index, N>{{ idx[n >= m_numIndices ? n : m_elements[which].representation[n]]... }}; 00094 } 00095 00096 template<typename Index> 00097 inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const 00098 { 00099 std::vector<Index> result; 00100 result.reserve(idx.size()); 00101 for (auto k : m_elements[which].representation) 00102 result.push_back(idx[k]); 00103 for (std::size_t i = m_numIndices; i < idx.size(); i++) 00104 result.push_back(idx[i]); 00105 return result; 00106 } 00107 00108 inline GroupElement ge(Generator const& g) const 00109 { 00110 GroupElement result; 00111 result.representation.reserve(m_numIndices); 00112 result.flags = g.flags; 00113 for (std::size_t k = 0; k < m_numIndices; k++) { 00114 if (k == (std::size_t)g.one) 00115 result.representation.push_back(g.two); 00116 else if (k == (std::size_t)g.two) 00117 result.representation.push_back(g.one); 00118 else 00119 result.representation.push_back(int(k)); 00120 } 00121 return result; 00122 } 00123 00124 GroupElement mul(GroupElement, GroupElement) const; 00125 inline GroupElement mul(Generator g1, GroupElement g2) const 00126 { 00127 return mul(ge(g1), g2); 00128 } 00129 00130 inline GroupElement mul(GroupElement g1, Generator g2) const 00131 { 00132 return mul(g1, ge(g2)); 00133 } 00134 00135 inline GroupElement mul(Generator g1, Generator g2) const 00136 { 00137 return mul(ge(g1), ge(g2)); 00138 } 00139 00140 inline int findElement(GroupElement e) const 00141 { 00142 for (auto ee : m_elements) { 00143 if (ee.representation == e.representation) 00144 return ee.flags ^ e.flags; 00145 } 00146 return -1; 00147 } 00148 00149 void updateGlobalFlags(int flagDiffOfSameGenerator); 00150 }; 00151 00152 // dynamic symmetry group that auto-adds the template parameters in the constructor 00153 template<typename... Gen> 00154 class DynamicSGroupFromTemplateArgs : public DynamicSGroup 00155 { 00156 public: 00157 inline DynamicSGroupFromTemplateArgs() : DynamicSGroup() 00158 { 00159 add_all(internal::type_list<Gen...>()); 00160 } 00161 inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) { } 00162 inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) { } 00163 inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) { DynamicSGroup::operator=(o); return *this; } 00164 inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) { DynamicSGroup::operator=(o); return *this; } 00165 00166 private: 00167 template<typename Gen1, typename... GenNext> 00168 inline void add_all(internal::type_list<Gen1, GenNext...>) 00169 { 00170 add(Gen1()); 00171 add_all(internal::type_list<GenNext...>()); 00172 } 00173 00174 inline void add_all(internal::type_list<>) 00175 { 00176 } 00177 }; 00178 00179 inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const 00180 { 00181 eigen_internal_assert(g1.representation.size() == m_numIndices); 00182 eigen_internal_assert(g2.representation.size() == m_numIndices); 00183 00184 GroupElement result; 00185 result.representation.reserve(m_numIndices); 00186 for (std::size_t i = 0; i < m_numIndices; i++) { 00187 int v = g2.representation[g1.representation[i]]; 00188 eigen_assert(v >= 0); 00189 result.representation.push_back(v); 00190 } 00191 result.flags = g1.flags ^ g2.flags; 00192 return result; 00193 } 00194 00195 inline void DynamicSGroup::add(int one, int two, int flags) 00196 { 00197 eigen_assert(one >= 0); 00198 eigen_assert(two >= 0); 00199 eigen_assert(one != two); 00200 00201 if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) { 00202 std::size_t newNumIndices = (one > two) ? one : two + 1; 00203 for (auto& gelem : m_elements) { 00204 gelem.representation.reserve(newNumIndices); 00205 for (std::size_t i = m_numIndices; i < newNumIndices; i++) 00206 gelem.representation.push_back(i); 00207 } 00208 m_numIndices = newNumIndices; 00209 } 00210 00211 Generator g{one, two, flags}; 00212 GroupElement e = ge(g); 00213 00214 /* special case for first generator */ 00215 if (m_elements.size() == 1) { 00216 while (!e.isId()) { 00217 m_elements.push_back(e); 00218 e = mul(e, g); 00219 } 00220 00221 if (e.flags > 0) 00222 updateGlobalFlags(e.flags); 00223 00224 // only add in case we didn't have identity 00225 if (m_elements.size() > 1) 00226 m_generators.push_back(g); 00227 return; 00228 } 00229 00230 int p = findElement(e); 00231 if (p >= 0) { 00232 updateGlobalFlags(p); 00233 return; 00234 } 00235 00236 std::size_t coset_order = m_elements.size(); 00237 m_elements.push_back(e); 00238 for (std::size_t i = 1; i < coset_order; i++) 00239 m_elements.push_back(mul(m_elements[i], e)); 00240 m_generators.push_back(g); 00241 00242 std::size_t coset_rep = coset_order; 00243 do { 00244 for (auto g : m_generators) { 00245 e = mul(m_elements[coset_rep], g); 00246 p = findElement(e); 00247 if (p < 0) { 00248 // element not yet in group 00249 m_elements.push_back(e); 00250 for (std::size_t i = 1; i < coset_order; i++) 00251 m_elements.push_back(mul(m_elements[i], e)); 00252 } else if (p > 0) { 00253 updateGlobalFlags(p); 00254 } 00255 } 00256 coset_rep += coset_order; 00257 } while (coset_rep < m_elements.size()); 00258 } 00259 00260 inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator) 00261 { 00262 switch (flagDiffOfSameGenerator) { 00263 case 0: 00264 default: 00265 // nothing happened 00266 break; 00267 case NegationFlag: 00268 // every element is it's own negative => whole tensor is zero 00269 m_globalFlags |= GlobalZeroFlag; 00270 break; 00271 case ConjugationFlag: 00272 // every element is it's own conjugate => whole tensor is real 00273 m_globalFlags |= GlobalRealFlag; 00274 break; 00275 case (NegationFlag | ConjugationFlag): 00276 // every element is it's own negative conjugate => whole tensor is imaginary 00277 m_globalFlags |= GlobalImagFlag; 00278 break; 00279 /* NOTE: 00280 * since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator 00281 * causes the tensor to be real and the next one to be imaginary, this will 00282 * trivially give the correct result 00283 */ 00284 } 00285 } 00286 00287 } // end namespace Eigen 00288 00289 #endif // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H 00290 00291 /* 00292 * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle; 00293 */