v0.7.26
Tensor2_symmetric_value.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 /* A general version, not for pointers. */
4 
5 #include <ostream>
6 #ifdef FTENSOR_DEBUG
7 #include <sstream>
8 #include <stdexcept>
9 #endif
10 
11 namespace FTensor
12 {
13  template <class T, int Tensor_Dim> class Tensor2_symmetric
14  {
15  T data[(Tensor_Dim * (Tensor_Dim + 1)) / 2];
16 
17  public:
18  template <class... U> Tensor2_symmetric(U... d) : data{d...}
19  {
20  static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
21  "Incorrect number of Arguments. Constructor should "
22  "initialize the entire Tensor");
23  }
24 
26 
27  /* There are two operator(int,int)'s, one for non-consts that lets you
28  change the value, and one for consts that doesn't. */
29 
30  T &operator()(const int N1, const int N2)
31  {
32 #ifdef FTENSOR_DEBUG
33  if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
34  {
35  std::stringstream s;
36  s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
37  << ">.operator(" << N1 << "," << N2 << ")" << std::endl;
38  throw std::out_of_range(s.str());
39  }
40 #endif
41  return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
42  : data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
43  }
44 
45  T operator()(const int N1, const int N2) const
46  {
47 #ifdef FTENSOR_DEBUG
48  if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
49  {
50  std::stringstream s;
51  s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
52  << ">.operator(" << N1 << "," << N2 << ") const" << std::endl;
53  throw std::out_of_range(s.str());
54  }
55 #endif
56  return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
57  : data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
58  }
59 
60  /* These operator()'s are the first part in constructing template
61  expressions. They can be used to slice off lower dimensional
62  parts. They are not entirely safe, since you can accidentaly use a
63  higher dimension than what is really allowed (like Dim=5). */
64 
65  /* This returns a Tensor2_Expr, since the indices are not really
66  symmetric anymore since they cover different dimensions. */
67 
68  template <char i, char j, int Dim0, int Dim1>
69  typename std::enable_if<
70  (Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
71  Tensor2_Expr<Tensor2_symmetric<T, Tensor_Dim>, T, Dim0, Dim1, i, j>>::type
72  operator()(const Index<i, Dim0>, const Index<j, Dim1>)
73  {
74  return Tensor2_Expr<Tensor2_symmetric<T, Tensor_Dim>, T, Dim0, Dim1, i,
75  j>(*this);
76  }
77 
78  template <char i, char j, int Dim0, int Dim1>
79  typename std::enable_if<(Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
81  T, Dim0, Dim1, i, j>>::type
82  operator()(const Index<i, Dim0>, const Index<j, Dim1>) const
83  {
85  Dim1, i, j>(*this);
86  }
87 
88  /* This returns a Tensor2_symmetric_Expr, since the indices are still
89  symmetric on the lower dimensions. */
90 
91  template <char i, char j, int Dim>
92  typename std::enable_if<
93  (Tensor_Dim >= Dim),
95  operator()(const Index<i, Dim>, const Index<j, Dim>)
96  {
98  i, j>(*this);
99  }
100 
101  template <char i, char j, int Dim>
102  typename std::enable_if<
103  (Tensor_Dim >= Dim),
105  j>>::type
106  operator()(const Index<i, Dim>, const Index<j, Dim>) const
107  {
109  Dim, i, j>(*this);
110  }
111 
112  /* This is for expressions where a number is used for one slot, and
113  an index for another, yielding a Tensor1_Expr. The non-const
114  versions don't actually create a Tensor2_number_rhs_[01] object.
115  They create a Tensor1_Expr directly, which provides the
116  appropriate indexing operators. The const versions do create a
117  Tensor2_number_[01]. */
118 
119  template <char i, int N, int Dim>
120  typename std::enable_if<
121  (Tensor_Dim >= Dim && Tensor_Dim > N),
123  T, Dim, i>>::type
124  operator()(const Index<i, Dim>, const Number<N> &)
125  {
126  using TensorExpr
129  }
130 
131  template <char i, int N, int Dim>
132  typename std::enable_if<
133  (Tensor_Dim >= Dim && Tensor_Dim > N),
135  T, Dim, i>>::type
136  operator()(const Index<i, Dim>, const Number<N> &) const
137  {
138  using TensorExpr
140  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
141  }
142 
143  template <char i, int N, int Dim>
144  typename std::enable_if<
145  (Tensor_Dim > N && Tensor_Dim >= Dim),
147  T, Dim, i>>::type
148  operator()(const Number<N> &, const Index<i, Dim>)
149  {
150  using TensorExpr
153  }
154 
155  template <char i, int N, int Dim>
156  typename std::enable_if<
157  (Tensor_Dim > N && Tensor_Dim >= Dim),
159  T, Dim, i>>::type
160  operator()(const Number<N> &, const Index<i, Dim>) const
161  {
162  using TensorExpr
164  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
165  }
166 
167  /* Specializations for using actual numbers instead of Number<> */
168 
169  template <char i, int Dim>
170  typename std::enable_if<
171  (Tensor_Dim >= Dim),
173  T, Dim, i>>::type
174  operator()(const Index<i, Dim>, const int N) const
175  {
176  using TensorExpr
178  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
179  }
180 
181  template <char i, int Dim>
182  typename std::enable_if<
183  (Tensor_Dim >= Dim),
185  T, Dim, i>>::type
186  operator()(const int N, const Index<i, Dim>) const
187  {
188  using TensorExpr
190  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
191  }
192 
193  /* These two operator()'s return the Tensor2 with internal
194  contractions, yielding a T. I have to specify one for both
195  const and non-const because otherwise the compiler will use the
196  operator() which gives a Tensor2_Expr<>. */
197 
198  template <char i, int Dim>
199  typename std::enable_if<(Tensor_Dim >= Dim), T>::type
200  operator()(const Index<i, Dim>, const Index<i, Dim>)
201  {
202  return internal_contract(Number<Dim>());
203  }
204 
205  template <char i, int Dim>
206  typename std::enable_if<(Tensor_Dim >= Dim), T>::type
207  operator()(const Index<i, Dim>, const Index<i, Dim>) const
208  {
209  return internal_contract(Number<Dim>());
210  }
211 
212  private:
213  template <int N> T internal_contract(const Number<N> &) const
214  {
215  return data[N - 1 + ((N - 1) * (2 * Tensor_Dim - N)) / 2]
217  }
218 
219  T internal_contract(const Number<1> &) const { return data[0]; }
220  };
221 }
222 
223 /// JSON compatible output. It only outputs unique elements, so a 3x3
224 /// symmetric matrix only outputs 6 elements.
225 
226 namespace FTensor
227 {
228  template <class T, int Tensor_Dim>
230  std::ostream &os, const FTensor::Tensor2_symmetric<T, Tensor_Dim> &t,
231  const int &i)
232  {
233  os << '[';
234  for(int j = i; j + 1 < Tensor_Dim; ++j)
235  {
236  os << t(i, j) << ',';
237  }
238  if(Tensor_Dim > 0)
239  {
240  os << t(i, Tensor_Dim - 1);
241  }
242  os << ']';
243  return os;
244  }
245 }
246 
247 template <class T, int Tensor_Dim>
248 std::ostream &operator<<(std::ostream &os,
250 {
251  os << '[';
252  for(int i = 0; i + 1 < Tensor_Dim; ++i)
253  {
255  os << ',';
256  }
257  if(Tensor_Dim > 0)
258  {
259  FTensor::Tensor2_symmetric_ostream_row(os, t, Tensor_Dim - 1);
260  }
261  os << ']';
262  return os;
263 }
264 
265 namespace FTensor
266 {
267  template <class T, int Tensor_Dim>
268  std::istream &
271  const int &i)
272  {
273  char c;
274  is >> c;
275  for(int j = i; j + 1 < Tensor_Dim; ++j)
276  {
277  is >> t(i, j) >> c;
278  }
279  if(Tensor_Dim > 0)
280  {
281  is >> t(i, Tensor_Dim - 1);
282  }
283  is >> c;
284  return is;
285  }
286 }
287 
288 template <class T, int Tensor_Dim>
289 std::istream &
291 {
292  char c;
293  is >> c;
294  for(int i = 0; i + 1 < Tensor_Dim; ++i)
295  {
297  is >> c;
298  }
299  if(Tensor_Dim > 0)
300  {
301  FTensor::Tensor2_symmetric_istream_row(is, t, Tensor_Dim - 1);
302  }
303  is >> c;
304  return is;
305 }
T & operator()(const int N1, const int N2)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T *> &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
std::ostream & Tensor2_symmetric_ostream_row(std::ostream &os, const FTensor::Tensor2_symmetric< T, Tensor_Dim > &t, const int &i)
T internal_contract(const Number< 1 > &) const
T data[(Tensor_Dim *(Tensor_Dim+1))/2]
Fully Antisymmetric Levi-Civita Tensor.
T internal_contract(const Number< N > &) const
std::istream & Tensor2_symmetric_istream_row(std::istream &is, FTensor::Tensor2_symmetric< T, Tensor_Dim > &t, const int &i)
T operator()(const int N1, const int N2) const
std::istream & operator>>(std::istream &is, FTensor::Tensor2_symmetric< T, Tensor_Dim > &t)
const int N
Definition: speed_test.cpp:3
std::ostream & operator<<(std::ostream &os, const FTensor::Tensor2_symmetric< T, Tensor_Dim > &t)