src/TensorSymmetry/Symmetry.h

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_SYMMETRY_H
0011 #define EIGEN_CXX11_TENSORSYMMETRY_SYMMETRY_H
0012
0013 namespace Eigen {
0014
0015 enum {
0016   NegationFlag           = 0x01,
0017   ConjugationFlag        = 0x02
0018 };
0019
0020 enum {
0021   GlobalRealFlag         = 0x01,
0022   GlobalImagFlag         = 0x02,
0023   GlobalZeroFlag         = 0x03
0024 };
0025
0026 namespace internal {
0027
0028 template<std::size_t NumIndices, typename... Sym>                   struct tensor_symmetry_pre_analysis;
0029 template<std::size_t NumIndices, typename... Sym>                   struct tensor_static_symgroup;
0030 template<bool instantiate, std::size_t NumIndices, typename... Sym> struct tensor_static_symgroup_if;
0031 template<typename Tensor_> struct tensor_symmetry_calculate_flags;
0032 template<typename Tensor_> struct tensor_symmetry_assign_value;
0033 template<typename... Sym> struct tensor_symmetry_num_indices;
0034
0035 } // end namespace internal
0036
0037 template<int One_, int Two_>
0038 struct Symmetry
0039 {
0040   static_assert(One_ != Two_, "Symmetries must cover distinct indices.");
0041   constexpr static int One = One_;
0042   constexpr static int Two = Two_;
0043   constexpr static int Flags = 0;
0044 };
0045
0046 template<int One_, int Two_>
0047 struct AntiSymmetry
0048 {
0049   static_assert(One_ != Two_, "Symmetries must cover distinct indices.");
0050   constexpr static int One = One_;
0051   constexpr static int Two = Two_;
0052   constexpr static int Flags = NegationFlag;
0053 };
0054
0055 template<int One_, int Two_>
0056 struct Hermiticity
0057 {
0058   static_assert(One_ != Two_, "Symmetries must cover distinct indices.");
0059   constexpr static int One = One_;
0060   constexpr static int Two = Two_;
0061   constexpr static int Flags = ConjugationFlag;
0062 };
0063
0064 template<int One_, int Two_>
0065 struct AntiHermiticity
0066 {
0067   static_assert(One_ != Two_, "Symmetries must cover distinct indices.");
0068   constexpr static int One = One_;
0069   constexpr static int Two = Two_;
0070   constexpr static int Flags = ConjugationFlag | NegationFlag;
0071 };
0072
0073 /** \class DynamicSGroup
0074   * \ingroup TensorSymmetry_Module
0075   *
0076   * \brief Dynamic symmetry group
0077   *
0078   * The %DynamicSGroup class represents a symmetry group that need not be known at
0079   * compile time. It is useful if one wants to support arbitrary run-time defineable
0080   * symmetries for tensors, but it is also instantiated if a symmetry group is defined
0081   * at compile time that would be either too large for the compiler to reasonably
0082   * generate (using templates to calculate this at compile time is very inefficient)
0083   * or that the compiler could generate the group but that it wouldn't make sense to
0084   * unroll the loop for setting coefficients anymore.
0085   */
0086 class DynamicSGroup;
0087
0088 /** \internal
0089   *
0090   * \class DynamicSGroupFromTemplateArgs
0091   * \ingroup TensorSymmetry_Module
0092   *
0093   * \brief Dynamic symmetry group, initialized from template arguments
0094   *
0095   * This class is a child class of DynamicSGroup. It uses the template arguments
0096   * specified to initialize itself.
0097   */
0098 template<typename... Gen>
0099 class DynamicSGroupFromTemplateArgs;
0100
0101 /** \class StaticSGroup
0102   * \ingroup TensorSymmetry_Module
0103   *
0104   * \brief Static symmetry group
0105   *
0106   * This class represents a symmetry group that is known and resolved completely
0107   * at compile time. Ideally, no run-time penalty is incurred compared to the
0108   * manual unrolling of the symmetry.
0109   *
0110   * <b><i>CAUTION:</i></b>
0111   *
0112   * Do not use this class directly for large symmetry groups. The compiler
0113   * may run into a limit, or segfault or in the very least will take a very,
0114   * very, very long time to compile the code. Use the SGroup class instead
0115   * if you want a static group. That class contains logic that will
0116   * automatically select the DynamicSGroup class instead if the symmetry
0117   * group becomes too large. (In that case, unrolling may not even be
0118   * beneficial.)
0119   */
0120 template<typename... Gen>
0121 class StaticSGroup;
0122
0123 /** \class SGroup
0124   * \ingroup TensorSymmetry_Module
0125   *
0126   * \brief Symmetry group, initialized from template arguments
0127   *
0128   * This class represents a symmetry group whose generators are already
0129   * known at compile time. It may or may not be resolved at compile time,
0130   * depending on the estimated size of the group.
0131   *
0132   * \sa StaticSGroup
0133   * \sa DynamicSGroup
0134   */
0135 template<typename... Gen>
0136 class SGroup : public internal::tensor_symmetry_pre_analysis<internal::tensor_symmetry_num_indices<Gen...>::value, Gen...>::root_type
0137 {
0138   public:
0139     constexpr static std::size_t NumIndices = internal::tensor_symmetry_num_indices<Gen...>::value;
0140     typedef typename internal::tensor_symmetry_pre_analysis<NumIndices, Gen...>::root_type Base;
0141
0142     // make standard constructors + assignment operators public
0143     inline SGroup() : Base() { }
0144     inline SGroup(const SGroup<Gen...>& other) : Base(other) { }
0145     inline SGroup(SGroup<Gen...>&& other) : Base(other) { }
0146     inline SGroup<Gen...>& operator=(const SGroup<Gen...>& other) { Base::operator=(other); return *this; }
0147     inline SGroup<Gen...>& operator=(SGroup<Gen...>&& other) { Base::operator=(other); return *this; }
0148
0149     // all else is defined in the base class
0150 };
0151
0152 namespace internal {
0153
0154 template<typename... Sym> struct tensor_symmetry_num_indices
0155 {
0156   constexpr static std::size_t value = 1;
0157 };
0158
0159 template<int One_, int Two_, typename... Sym> struct tensor_symmetry_num_indices<Symmetry<One_, Two_>, Sym...>
0160 {
0161 private:
0162   constexpr static std::size_t One = static_cast<std::size_t>(One_);
0163   constexpr static std::size_t Two = static_cast<std::size_t>(Two_);
0164   constexpr static std::size_t Three = tensor_symmetry_num_indices<Sym...>::value;
0165
0166   // don't use std::max, since it's not constexpr until C++14...
0167   constexpr static std::size_t maxOneTwoPlusOne = ((One > Two) ? One : Two) + 1;
0168 public:
0169   constexpr static std::size_t value = (maxOneTwoPlusOne > Three) ? maxOneTwoPlusOne : Three;
0170 };
0171
0172 template<int One_, int Two_, typename... Sym> struct tensor_symmetry_num_indices<AntiSymmetry<One_, Two_>, Sym...>
0173   : public tensor_symmetry_num_indices<Symmetry<One_, Two_>, Sym...> {};
0174 template<int One_, int Two_, typename... Sym> struct tensor_symmetry_num_indices<Hermiticity<One_, Two_>, Sym...>
0175   : public tensor_symmetry_num_indices<Symmetry<One_, Two_>, Sym...> {};
0176 template<int One_, int Two_, typename... Sym> struct tensor_symmetry_num_indices<AntiHermiticity<One_, Two_>, Sym...>
0177   : public tensor_symmetry_num_indices<Symmetry<One_, Two_>, Sym...> {};
0178
0179 /** \internal
0180   *
0181   * \class tensor_symmetry_pre_analysis
0182   * \ingroup TensorSymmetry_Module
0183   *
0184   * \brief Pre-select whether to use a static or dynamic symmetry group
0185   *
0186   * When a symmetry group could in principle be determined at compile time,
0187   * this template implements the logic whether to actually do that or whether
0188   * to rather defer that to runtime.
0189   *
0190   * The logic is as follows:
0191   * <dl>
0192   * <dt><b>No generators (trivial symmetry):</b></dt>
0193   * <dd>Use a trivial static group. Ideally, this has no performance impact
0194   *     compared to not using symmetry at all. In practice, this might not
0195   *     be the case.</dd>
0196   * <dt><b>More than 4 generators:</b></dt>
0197   * <dd>Calculate the group at run time, it is likely far too large for the
0198   *     compiler to be able to properly generate it in a realistic time.</dd>
0199   * <dt><b>Up to and including 4 generators:</b></dt>
0200   * <dd>Actually enumerate all group elements, but then check how many there
0201   *     are. If there are more than 16, it is unlikely that unrolling the
0202   *     loop (as is done in the static compile-time case) is sensible, so
0203   *     use a dynamic group instead. If there are at most 16 elements, actually
0204   *     use that static group. Note that the largest group with 4 generators
0205   *     still compiles with reasonable resources.</dd>
0206   * </dl>
0207   *
0208   * Note: Example compile time performance with g++-4.6 on an Intenl Core i5-3470
0209   *       with 16 GiB RAM (all generators non-redundant and the subgroups don't
0210   *       factorize):
0211   *
0212   *          # Generators          -O0 -ggdb               -O2
0213   *          -------------------------------------------------------------------
0214   *          1                 0.5 s  /   250 MiB     0.45s /   230 MiB
0215   *          2                 0.5 s  /   260 MiB     0.5 s /   250 MiB
0216   *          3                 0.65s  /   310 MiB     0.62s /   310 MiB
0217   *          4                 2.2 s  /   860 MiB     1.7 s /   770 MiB
0218   *          5               130   s  / 13000 MiB   120   s / 11000 MiB
0219   *
0220   * It is clear that everything is still very efficient up to 4 generators, then
0221   * the memory and CPU requirements become unreasonable. Thus we only instantiate
0222   * the template group theory logic if the number of generators supplied is 4 or
0223   * lower, otherwise this will be forced to be done during runtime, where the
0224   * algorithm is reasonably fast.
0225   */
0226 template<std::size_t NumIndices>
0227 struct tensor_symmetry_pre_analysis<NumIndices>
0228 {
0229   typedef StaticSGroup<> root_type;
0230 };
0231
0232 template<std::size_t NumIndices, typename Gen_, typename... Gens_>
0233 struct tensor_symmetry_pre_analysis<NumIndices, Gen_, Gens_...>
0234 {
0235   constexpr static std::size_t max_static_generators = 4;
0236   constexpr static std::size_t max_static_elements = 16;
0237   typedef tensor_static_symgroup_if<(sizeof...(Gens_) + 1 <= max_static_generators), NumIndices, Gen_, Gens_...> helper;
0238   constexpr static std::size_t possible_size = helper::size;
0239
0240   typedef typename conditional<
0241     possible_size == 0 || possible_size >= max_static_elements,
0242     DynamicSGroupFromTemplateArgs<Gen_, Gens_...>,
0243     typename helper::type
0244   >::type root_type;
0245 };
0246
0247 template<bool instantiate, std::size_t NumIndices, typename... Gens>
0248 struct tensor_static_symgroup_if
0249 {
0250   constexpr static std::size_t size = 0;
0251   typedef void type;
0252 };
0253
0254 template<std::size_t NumIndices, typename... Gens>
0255 struct tensor_static_symgroup_if<true, NumIndices, Gens...> : tensor_static_symgroup<NumIndices, Gens...> {};
0256
0257 template<typename Tensor_>
0258 struct tensor_symmetry_assign_value
0259 {
0260   typedef typename Tensor_::Index Index;
0261   typedef typename Tensor_::Scalar Scalar;
0262   constexpr static std::size_t NumIndices = Tensor_::NumIndices;
0263
0264   static inline int run(const std::array<Index, NumIndices>& transformed_indices, int transformation_flags, int dummy, Tensor_& tensor, const Scalar& value_)
0265   {
0266     Scalar value(value_);
0267     if (transformation_flags & ConjugationFlag)
0268       value = numext::conj(value);
0269     if (transformation_flags & NegationFlag)
0270       value = -value;
0271     tensor.coeffRef(transformed_indices) = value;
0272     return dummy;
0273   }
0274 };
0275
0276 template<typename Tensor_>
0277 struct tensor_symmetry_calculate_flags
0278 {
0279   typedef typename Tensor_::Index Index;
0280   constexpr static std::size_t NumIndices = Tensor_::NumIndices;
0281
0282   static inline int run(const std::array<Index, NumIndices>& transformed_indices, int transform_flags, int current_flags, const std::array<Index, NumIndices>& orig_indices)
0283   {
0284     if (transformed_indices == orig_indices) {
0285       if (transform_flags & (ConjugationFlag | NegationFlag))
0286         return current_flags | GlobalImagFlag; // anti-hermitian diagonal
0287       else if (transform_flags & ConjugationFlag)
0288         return current_flags | GlobalRealFlag; // hermitian diagonal
0289       else if (transform_flags & NegationFlag)
0290         return current_flags | GlobalZeroFlag; // anti-symmetric diagonal
0291     }
0292     return current_flags;
0293   }
0294 };
0295
0296 template<typename Tensor_, typename Symmetry_, int Flags = 0>
0297 class tensor_symmetry_value_setter
0298 {
0299   public:
0300     typedef typename Tensor_::Index Index;
0301     typedef typename Tensor_::Scalar Scalar;
0302     constexpr static std::size_t NumIndices = Tensor_::NumIndices;
0303
0304     inline tensor_symmetry_value_setter(Tensor_& tensor, Symmetry_ const& symmetry, std::array<Index, NumIndices> const& indices)
0305       : m_tensor(tensor), m_symmetry(symmetry), m_indices(indices) { }
0306
0307     inline tensor_symmetry_value_setter<Tensor_, Symmetry_, Flags>& operator=(Scalar const& value)
0308     {
0309       doAssign(value);
0310       return *this;
0311     }
0312   private:
0313     Tensor_& m_tensor;
0314     Symmetry_ m_symmetry;
0315     std::array<Index, NumIndices> m_indices;
0316
0317     inline void doAssign(Scalar const& value)
0318     {
0319       #ifdef EIGEN_TENSOR_SYMMETRY_CHECK_VALUES
0320         int value_flags = m_symmetry.template apply<internal::tensor_symmetry_calculate_flags<Tensor_>, int>(m_indices, m_symmetry.globalFlags(), m_indices);
0321         if (value_flags & GlobalRealFlag)
0322           eigen_assert(numext::imag(value) == 0);
0323         if (value_flags & GlobalImagFlag)
0324           eigen_assert(numext::real(value) == 0);
0325       #endif
0326       m_symmetry.template apply<internal::tensor_symmetry_assign_value<Tensor_>, int>(m_indices, 0, m_tensor, value);
0327     }
0328 };
0329
0330 } // end namespace internal
0331
0332 } // end namespace Eigen
0333
0334 #endif // EIGEN_CXX11_TENSORSYMMETRY_SYMMETRY_H
0335
0336 /*
0337  * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
0338  */