1 /*
2  *            Copyright 2009-2020 The VOTCA Development Team
3  *                       (http://www.votca.org)
4  *
5  *      Licensed under the Apache License, Version 2.0 (the "License")
6  *
7  * You may not use this file except in compliance with the License.
8  * You may obtain a copy of the License at
9  *
10  *              http://www.apache.org/licenses/LICENSE-2.0
11  *
12  * Unless required by applicable law or agreed to in writing, software
13  * distributed under the License is distributed on an "AS IS" BASIS,
14  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  * See the License for the specific language governing permissions and
16  * limitations under the License.
17  *
18  */
19 
20 #ifndef VOTCA_TOOLS_NDIMVECTOR_H
21 #define VOTCA_TOOLS_NDIMVECTOR_H
22 
23 // Standard includes
24 #include "types.h"
25 #include <algorithm>
26 #include <array>
27 #include <cassert>
28 #include <numeric>
29 #include <vector>
30 
31 namespace votca {
32 namespace tools {
33 
34 /**
35  * \brief N-Dim Vector
36  *
37  * This is basically a normal std::vector with an additional
38  * operator(i,j,k,l...), which allows you to treat it as a higher dim object
39  * It is column major (the first index is the fast one)
40  *
41  */
42 
43 template <class T, int dim>
44 class NDimVector {
45 
46   static constexpr int dim_ = dim;
47 
48   template <typename... IndexTypes>
linearIndex(Index firstDimension,IndexTypes...otherDimensions)49   Index linearIndex(Index firstDimension, IndexTypes... otherDimensions) const {
50     static_assert((sizeof...(otherDimensions) + 1 == dim_),
51                   "Number of dimensions given does not match rank");
52     std::array<Index, dim_> indices = {firstDimension, otherDimensions...};
53     assert(std::all_of(indices.begin(), indices.end(),
54                        [](Index size) { return size >= 0; }) &&
55            "All indeces must be non-negative");
56     assert(std::equal(indices.begin(), indices.end(), dimensions_.begin(),
57                       [](Index index, Index dimension) {
58                         return index < dimension;
59                       }) &&
60            "All indeces must be smaller than dimension");
61     return std::inner_product(indices.begin(), indices.end(), offsets_.begin(),
62                               0);
63   }
64 
65  public:
rank()66   constexpr Index rank() { return dimensions_.size(); }
67 
68   NDimVector() = default;
69 
70   template <typename... IndexTypes>
NDimVector(Index firstDimension,IndexTypes...otherDimensions)71   NDimVector(Index firstDimension, IndexTypes... otherDimensions) {
72     // The number of dimensions used to construct a tensor must be equal to the
73     // rank of the tensor.
74     static_assert((sizeof...(otherDimensions) + 1 == dim_),
75                   "Number of dimensions given does not match rank");
76     dimensions_ = {firstDimension, otherDimensions...};
77 
78     assert(std::all_of(dimensions_.begin(), dimensions_.end(),
79                        [](Index size) { return size >= 0; }) &&
80            "All dimensions must be larger 0");
81     Index size = std::accumulate(dimensions_.begin(), dimensions_.end(), 1,
82                                  std::multiplies<Index>());
83     storage_.resize(size);
84     offsets_[0] = 1;
85     std::partial_sum(dimensions_.begin(), dimensions_.end() - 1,
86                      offsets_.begin() + 1, std::multiplies<Index>());
87   }
88 
89   template <typename... IndexTypes>
operator()90   T& operator()(Index firstDimension, IndexTypes... otherDimensions) {
91     return storage_[linearIndex(firstDimension, otherDimensions...)];
92     static_assert((sizeof...(otherDimensions) + 1 == dim_),
93                   "Number of dimensions given does not match rank");
94     std::array<Index, dim_> indices = {firstDimension, otherDimensions...};
95     Index linear_index =
96         std::inner_product(indices.begin(), indices.end(), offsets_.begin(), 0);
97     return storage_[linear_index];
98   }
99 
100   template <typename... IndexTypes>
operator()101   const T& operator()(Index firstDimension,
102                       IndexTypes... otherDimensions) const {
103     return storage_[linearIndex(firstDimension, otherDimensions...)];
104   }
105 
dimension(Index i)106   Index dimension(Index i) const { return dimensions_[i]; }
size()107   Index size() const { return storage_.size(); }
108 
begin()109   typename std::vector<T>::iterator begin() { return storage_.begin(); }
end()110   typename std::vector<T>::iterator end() { return storage_.end(); }
begin()111   typename std::vector<T>::const_iterator begin() const {
112     return storage_.begin();
113   }
end()114   typename std::vector<T>::const_iterator end() const { return storage_.end(); }
115 
116  private:
117   std::vector<T> storage_;
118   std::array<Index, dim> dimensions_;
119   std::array<Index, dim> offsets_;
120 };
121 }  // namespace tools
122 }  // namespace votca
123 
124 #endif  // VOTCA_TOOLS_NDIMVECTOR_H
125