v0.7.2
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>
15  {
16  T data[(Tensor_Dim*(Tensor_Dim+1))/2];
17  public:
19 
20  /* Tensor_Dim=1 */
22  {
24  }
25 
26  /* Tensor_Dim=2 */
27  Tensor2_symmetric(T d00, T d01, T d11)
28  {
30  }
31 
32  /* Tensor_Dim=3 */
33  Tensor2_symmetric(T d00, T d01, T d02, T d11, T d12, T d22)
34  {
35  Tensor2_symmetric_constructor<T,Tensor_Dim>(data,d00,d01,d02,d11,d12,d22);
36  }
37 
38  /* Tensor_Dim=4 */
39  Tensor2_symmetric(T d00, T d01, T d02, T d03, T d11, T d12, T d13, T d22,
40  T d23, T d33)
41  {
43  (data,d00,d01,d02,d03,d11,d12,d13,d22,d23,d33);
44  }
45 
46  /* There are two operator(int,int)'s, one for non-consts that lets you
47  change the value, and one for consts that doesn't. */
48 
49  T & operator()(const int N1, const int N2)
50  {
51 #ifdef FTENSOR_DEBUG
52  if(N1>=Tensor_Dim || N1<0 || N2>=Tensor_Dim || N2<0)
53  {
54  std::stringstream s;
55  s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
56  << ">.operator(" << N1 << "," << N2 << ")"
57  << std::endl;
58  throw std::runtime_error(s.str());
59  }
60 #endif
61  return N1>N2 ? data[N1+(N2*(2*Tensor_Dim-N2-1))/2]
62  : data[N2+(N1*(2*Tensor_Dim-N1-1))/2];
63  }
64 
65  T operator()(const int N1, const int N2) const
66  {
67 #ifdef FTENSOR_DEBUG
68  if(N1>=Tensor_Dim || N1<0 || N2>=Tensor_Dim || N2<0)
69  {
70  std::stringstream s;
71  s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
72  << ">.operator(" << N1 << "," << N2 << ") const"
73  << std::endl;
74  throw std::runtime_error(s.str());
75  }
76 #endif
77  return N1>N2 ? data[N1+(N2*(2*Tensor_Dim-N2-1))/2]
78  : data[N2+(N1*(2*Tensor_Dim-N1-1))/2];
79  }
80 
81  /* These operator()'s are the first part in constructing template
82  expressions. They can be used to slice off lower dimensional
83  parts. They are not entirely safe, since you can accidently use a
84  higher dimension than what is really allowed (like Dim=5). */
85 
86  /* This returns a Tensor2_Expr, since the indices are not really
87  symmetric anymore since they cover different dimensions. */
88 
89  template<char i, char j, int Dim0, int Dim1>
90  typename std::enable_if<(Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
91  Tensor2_Expr<Tensor2_symmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j> >::type
92  operator()(const Index<i,Dim0> , const Index<j,Dim1> )
93  {
94  return Tensor2_Expr<Tensor2_symmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j>
95  (*this);
96  }
97 
98  template<char i, char j, int Dim0, int Dim1>
99  typename std::enable_if<(Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
100  Tensor2_Expr<const Tensor2_symmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j> >::type
101  operator()(const Index<i,Dim0> , const Index<j,Dim1> ) const
102  {
104  (*this);
105  }
106 
107  /* This returns a Tensor2_symmetric_Expr, since the indices are still
108  symmetric on the lower dimensions. */
109 
110  template<char i, char j, int Dim>
111  typename std::enable_if<(Tensor_Dim >= Dim),
113  operator()(const Index<i,Dim> , const Index<j,Dim> )
114  {
116  (*this);
117  }
118 
119  template<char i, char j, int Dim>
120  typename std::enable_if<(Tensor_Dim >= Dim),
122  operator()(const Index<i,Dim> , const Index<j,Dim> ) const
123  {
125  T,Dim,i,j>(*this);
126  }
127 
128  /* This is for expressions where a number is used for one slot, and
129  an index for another, yielding a Tensor1_Expr. The non-const
130  versions don't actually create a Tensor2_number_rhs_[01] object.
131  They create a Tensor1_Expr directly, which provides the
132  appropriate indexing operators. The const versions do create a
133  Tensor2_number_[01]. */
134 
135  template<char i, int N, int Dim>
136  typename std::enable_if<(Tensor_Dim >= Dim && Tensor_Dim > N),
138  T,Dim,i> >::type
139  operator()(const Index<i,Dim> , const Number<N> &)
140  {
142  TensorExpr;
143  return Tensor1_Expr<TensorExpr,T,Dim,i>(*this);
144  }
145 
146  template<char i, int N, int Dim>
147  typename std::enable_if<(Tensor_Dim >= Dim && Tensor_Dim > N),
149  T,N>,T,Dim,i> >::type
150  operator()(const Index<i,Dim> , const Number<N> &) const
151  {
153  TensorExpr;
154  return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this));
155  }
156 
157  template<char i, int N, int Dim>
158  typename std::enable_if<(Tensor_Dim > N && Tensor_Dim >= Dim),
160  T,Dim,i> >::type
161  operator()(const Number<N> &, const Index<i,Dim> )
162  {
164  TensorExpr;
165  return Tensor1_Expr<TensorExpr,T,Dim,i>(*this);
166  }
167 
168  template<char i, int N, int Dim>
169  typename std::enable_if<(Tensor_Dim > N && Tensor_Dim >= Dim),
171  T,N>,T,Dim,i> >::type
172  operator()(const Number<N> &, const Index<i,Dim> ) const
173  {
175  TensorExpr;
176  return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this));
177  }
178 
179  /* Specializations for using actual numbers instead of Number<> */
180 
181  template<char i, int Dim>
182  typename std::enable_if<(Tensor_Dim >= Dim),
184  T>,T,Dim,i> >::type
185  operator()(const Index<i,Dim> , const int N) const
186  {
188  TensorExpr;
189  return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this,N));
190  }
191 
192  template<char i, int Dim>
193  typename std::enable_if<(Tensor_Dim >= Dim),
195  T>,T,Dim,i> >::type
196  operator()(const int N, const Index<i,Dim> ) const
197  {
199  TensorExpr;
200  return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this,N));
201  }
202 
203  /* These two operator()'s return the Tensor2 with internal
204  contractions, yielding a T. I have to specify one for both
205  const and non-const because otherwise the compiler will use the
206  operator() which gives a Tensor2_Expr<>. */
207 
208  template<char i, int Dim>
209  typename std::enable_if<(Tensor_Dim >= Dim),
210  T>::type operator()(const Index<i,Dim> , const Index<i,Dim> )
211  {
212  return internal_contract(Number<Dim>());
213  }
214 
215  template<char i, int Dim>
216  typename std::enable_if<(Tensor_Dim >= Dim),
217  T>::type operator()(const Index<i,Dim> , const Index<i,Dim> ) const
218  {
219  return internal_contract(Number<Dim>());
220  }
221  private:
222  template<int N>
223  T internal_contract(const Number<N> &) const
224  {
225  return data[N-1+((N-1)*(2*Tensor_Dim-N))/2]
227  }
228 
229  T internal_contract(const Number<1> &) const
230  {
231  return data[0];
232  }
233  };
234 }
235 
236 /// JSON compatible output. It only outputs unique elments, so a 3x3
237 /// symmetric matrix only outputs 6 elements.
238 
239 namespace FTensor
240 {
241  template <class T, int Tensor_Dim>
242  std::ostream &
244  (std::ostream &os,
246  const int &i)
247  {
248  os << '[';
249  for(int j=i; j+1<Tensor_Dim; ++j)
250  { os << t(i,j) << ','; }
251  if(Tensor_Dim>0)
252  { os << t(i,Tensor_Dim-1); }
253  os << ']';
254  return os;
255  }
256 }
257 
258 template <class T, int Tensor_Dim>
259 std::ostream & operator<<(std::ostream &os,
261 {
262  os << '[';
263  for(int i=0; i+1<Tensor_Dim; ++i)
264  {
266  os << ',';
267  }
268  if(Tensor_Dim>0)
269  {
270  FTensor::Tensor2_symmetric_ostream_row(os,t,Tensor_Dim-1);
271  }
272  os << ']';
273  return os;
274 }
275 
276 
277 namespace FTensor
278 {
279  template <class T, int Tensor_Dim>
280  std::istream &
282  (std::istream &is,
284  const int &i)
285  {
286  char c;
287  is >> c;
288  for(int j=i; j+1<Tensor_Dim; ++j)
289  { is >> t(i,j) >> c; }
290  if(Tensor_Dim>0)
291  { is >> t(i,Tensor_Dim-1); }
292  is >> c;
293  return is;
294  }
295 }
296 
297 template <class T, int Tensor_Dim>
298 std::istream & operator>>(std::istream &is,
300 {
301  char c;
302  is >> c;
303  for(int i=0; i+1<Tensor_Dim; ++i)
304  {
306  is >> c;
307  }
308  if(Tensor_Dim>0)
309  {
310  FTensor::Tensor2_symmetric_istream_row(is,t,Tensor_Dim-1);
311  }
312  is >> c;
313  return is;
314 }
T & operator()(const int N1, const int N2)
std::ostream & Tensor2_symmetric_ostream_row(std::ostream &os, const FTensor::Tensor2_symmetric< T, Tensor_Dim > &t, const int &i)
Tensor2_symmetric(T d00, T d01, T d02, T d03, T d11, T d12, T d13, T d22, T d23, T d33)
T internal_contract(const Number< 1 > &) const
T data[(Tensor_Dim *(Tensor_Dim+1))/2]
JSON compatible output.
Tensor2_symmetric(T d00, T d01, T d02, T d11, T d12, T d22)
T internal_contract(const Number< N > &) const
Tensor2_symmetric(T d00, T d01, T d11)
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)