v0.9.0
Tensor3_value.hpp
Go to the documentation of this file.
1 /* A general version, not for pointers. */
2 
3 #pragma once
4 
5 #include "Tensor3_contracted.hpp"
6 
7 #include <ostream>
8 
9 namespace FTensor
10 {
11  template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
12  class Tensor3
13  {
14  T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2];
15 
16  public:
17  template <class... U> Tensor3(U... d) : data{d...}
18  {
19  static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
20  "Incorrect number of Arguments. Constructor should "
21  "initialize the entire Tensor");
22  }
23 
24  Tensor3() {}
25 
26  /* There are two operator(int,int,int)'s, one for non-consts that lets you
27  change the value, and one for consts that doesn't. */
28 
29  T &operator()(const int N1, const int N2, const int N3)
30  {
31 #ifdef FTENSOR_DEBUG
32  if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
33  || N3 >= Tensor_Dim2 || N3 < 0)
34  {
35  std::stringstream s;
36  s << "Bad index in Tensor3<T," << Tensor_Dim0 << "," << Tensor_Dim1
37  << "," << Tensor_Dim2 << ">.operator(" << N1 << "," << N2 << ","
38  << N3 << ")" << std::endl;
39  throw std::out_of_range(s.str());
40  }
41 #endif
42  return data[N1][N2][N3];
43  }
44 
45  T operator()(const int N1, const int N2, const int N3) const
46  {
47 #ifdef FTENSOR_DEBUG
48  if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
49  || N3 >= Tensor_Dim2 || N3 < 0)
50  {
51  std::stringstream s;
52  s << "Bad index in Tensor3<T," << Tensor_Dim0 << "," << Tensor_Dim1
53  << "," << Tensor_Dim2 << ">.operator(" << N1 << "," << N2 << ","
54  << N3 << ") const" << std::endl;
55  throw std::out_of_range(s.str());
56  }
57 #endif
58  return data[N1][N2][N3];
59  }
60 
61  /* These operator()'s are the first part in constructing template
62  expressions. */
63 
64  template <char i, char j, char k, int Dim0, int Dim1, int Dim2>
65  typename std::enable_if<
66  (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
68  Dim1, Dim2, i, j, k>>::type
69  operator()(const Index<i, Dim0>, const Index<j, Dim1>,
70  const Index<k, Dim2>)
71  {
73  Dim0, Dim1, Dim2, i, j, k>(*this);
74  }
75 
76  template <char i, char j, char k, int Dim0, int Dim1, int Dim2>
77  typename std::enable_if<
78  (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
80  Dim0, Dim1, Dim2, i, j, k>>::type
81  operator()(const Index<i, Dim0>, const Index<j, Dim1>,
82  const Index<k, Dim2>) const
83  {
84  return Tensor3_Expr<
86  Dim2, i, j, k>(*this);
87  }
88 
89  /* These operators are for internal contractions. We can not do
90  something like A(i,j,j) where i and j have different dimensions,
91  because it becomes ambiguous. */
92  /*TODO I dont know how i and j having different dimensions can be
93  * ambiguous?? Since you do Use that in A(j,j,i)*/
94  /* A(i,j,j) */
95 
96  template <char i, char j, int Dim>
97  typename std::enable_if<
98  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim),
102  T, Dim, i>>::type
103  operator()(const Index<i, Dim>, const Index<j, Dim>,
104  const Index<j, Dim>) const
105  {
106  using TensorExpr = Tensor3_contracted_12<
108  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
109  }
110 
111  /* A(j,i,j) */
112 
113  template <char i, char j, int Dim>
114  typename std::enable_if<
115  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim),
116  Tensor1_Expr<
119  T, Dim, i>>::type
120  operator()(const Index<j, Dim>, const Index<i, Dim>,
121  const Index<j, Dim>) const
122  {
123  using TensorExpr = Tensor3_contracted_02<
125  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
126  }
127 
128  /* A(j,j,i) */
129 
130  template <char i, char j, int Dim, int Dim2>
131  typename std::enable_if<
132  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim2),
133  Tensor1_Expr<
136  T, Dim2, i>>::type
137  operator()(const Index<j, Dim>, const Index<j, Dim>,
138  const Index<i, Dim2>) const
139  {
140  using TensorExpr = Tensor3_contracted_01<
142  return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
143  }
144 
145  /* This is for expressions where a number is used for one slot, and
146  indices for the others, yielding a Tensor2_Expr or
147  Tensor2_symmetric_Expr. The non-const versions don't actually
148  create a Tensor3_number_rhs_* object, while the const versions
149  do create a Tensor3_number_*. */
150 
151  /* First slot. */
152 
153  template <char i, char j, int N, int Dim1, int Dim2>
154  typename std::enable_if<
155  (Tensor_Dim0 > N && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
158  T, Dim1, Dim2, i, j>>::type
159  operator()(const Number<N>, const Index<i, Dim1>, const Index<j, Dim2>)
160  {
161  using TensorExpr = Tensor3_number_rhs_0<
164  }
165 
166  template <char i, char j, int N, int Dim1, int Dim2>
167  typename std::enable_if<
168  (Tensor_Dim0 > N && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
169  Tensor2_Expr<
172  T, Dim1, Dim2, i, j>>::type
173  operator()(const Number<N>, const Index<i, Dim1>,
174  const Index<j, Dim2>) const
175  {
176  using TensorExpr = Tensor3_number_0<
178  return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(TensorExpr(*this));
179  }
180 
181  /* Second slot. */
182 
183  template <char i, char j, int N, int Dim0, int Dim2>
184  typename std::enable_if<
185  (Tensor_Dim0 >= Dim0 && Tensor_Dim1 > N && Tensor_Dim2 >= Dim2),
188  T, Dim0, Dim2, i, j>>::type
189  operator()(const Index<i, Dim0>, const Number<N>, const Index<j, Dim2>)
190  {
191  using TensorExpr = Tensor3_number_rhs_0<
194  }
195 
196  template <char i, char j, int N, int Dim0, int Dim2>
197  typename std::enable_if<
198  (Tensor_Dim0 >= Dim0 && Tensor_Dim1 > N && Tensor_Dim2 >= Dim2),
199  Tensor2_Expr<
202  T, Dim0, Dim2, i, j>>::type
203  operator()(const Index<i, Dim0>, const Number<N>,
204  const Index<j, Dim2>) const
205  {
206  using TensorExpr = Tensor3_number_0<
208  return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(TensorExpr(*this));
209  }
210 
211  /* Third slot. */
212  // TODO allow two different dimensions here
213  template <char i, char j, int N, int Dim>
214  typename std::enable_if<
215  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 > N),
218  T, N>,
219  T, Dim, i, j>>::type
220  operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>)
221  {
222  using TensorExpr = Tensor3_number_rhs_2<
225  }
226  // TODO allow two different dimensions here
227  template <char i, char j, int N, int Dim>
228  typename std::enable_if<
229  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 > N),
233  T, Dim, i, j>>::type
234  operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>) const
235  {
236  using TensorExpr = Tensor3_number_2<
239  TensorExpr(*this));
240  }
241 
242  /* This is for expressions where a number is used for two slots, and
243  an Index for the other, yielding a Tensor1_Expr. The non-const
244  versions don't actually create a Tensor3_number_rhs_* object,
245  while the const versions do create a Tensor3_number_*. */
246 
247  /* Index in first slot. */
248 
249  template <char i, int N1, int N2, int Dim>
250  typename std::enable_if<
251  (Tensor_Dim0 >= Dim && Tensor_Dim1 > N1 && Tensor_Dim2 > N2),
252  Tensor1_Expr<
255  T, Dim, i>>::type
256  operator()(const Index<i, Dim> index, const Number<N1>, const Number<N2>)
257  {
258  using TensorExpr = Tensor3_number_rhs_12<
261  }
262 
263  template <char i, int N1, int N2, int Dim>
264  typename std::enable_if<
265  (Tensor_Dim0 >= Dim && Tensor_Dim1 > N1 && Tensor_Dim2 > N2),
266  Tensor1_Expr<
269  T, Dim, i>>::type
270  operator()(const Index<i, Dim> index, const Number<N1>,
271  const Number<N2>) const
272  {
273  using TensorExpr = Tensor3_number_12<
275  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
276  }
277 
278  /* Index in second slot. I use the same structures as for the Index
279  in the first slot since the tensor is symmetric on the first two
280  indices. */
281 
282  template <char i, int N0, int N2, int Dim>
283  typename std::enable_if<
284  (Tensor_Dim0 > N0 && Tensor_Dim1 >= Dim && Tensor_Dim2 > N2),
285  Tensor1_Expr<
288  T, Dim, i>>::type
289  operator()(const Number<N0>, const Index<i, Dim> index, const Number<N2>)
290  {
291  using TensorExpr = Tensor3_number_rhs_12<
294  }
295 
296  template <char i, int N0, int N2, int Dim>
297  typename std::enable_if<
298  (Tensor_Dim0 > N0 && Tensor_Dim1 >= Dim && Tensor_Dim2 > N2),
299  Tensor1_Expr<
302  T, Dim, i>>::type
303  operator()(const Number<N0>, const Index<i, Dim> index,
304  const Number<N2>) const
305  {
306  using TensorExpr = Tensor3_number_12<
308  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
309  }
310 
311  /* Index in third slot. */
312 
313  template <char i, int N0, int N1, int Dim>
314  typename std::enable_if<
315  (Tensor_Dim0 > N0 && Tensor_Dim1 > N1 && Tensor_Dim2 >= Dim),
316  Tensor1_Expr<
319  T, Dim, i>>::type
320  operator()(const Number<N0>, const Number<N1>, const Index<i, Dim> index)
321  {
322  using TensorExpr = Tensor3_number_rhs_01<
325  }
326 
327  template <char i, int N0, int N1, int Dim>
328  typename std::enable_if<
329  (Tensor_Dim0 > N0 && Tensor_Dim1 > N1 && Tensor_Dim2 >= Dim),
330  Tensor1_Expr<
333  T, Dim, i>>::type
334  operator()(const Number<N0>, const Number<N1>,
335  const Index<i, Dim> index) const
336  {
337  using TensorExpr = Tensor3_number_01<
339  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
340  }
341 
342  /* Specializations for using actual numbers instead of Number<>.
343  This is for expressions where an actual number is used for one
344  slot, and indices for the others, yielding a Tensor2_Expr or
345  Tensor2_symmetric_Expr. I only define the const versions. */
346 
347  /* First slot. */
348 
349  template <char i, char j, int Dim1, int Dim2>
350  typename std::enable_if<
351  (Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
354  T, Dim1, Dim2, i, j>>::type
355  operator()(const int N, const Index<i, Dim1>, const Index<j, Dim2>) const
356  {
357  using TensorExpr = Tensor3_numeral_0<
360  TensorExpr(*this, N));
361  }
362 
363  /* Second slot. */
364 
365  template <char i, char j, int Dim0, int Dim2>
366  typename std::enable_if<
367  (Tensor_Dim0 >= Dim0 && Tensor_Dim2 >= Dim2),
370  T, Dim0, Dim2, i, j>>::type
371  operator()(const Index<i, Dim0>, const int N, const Index<j, Dim2>) const
372  {
373  using TensorExpr = Tensor3_numeral_0<
376  TensorExpr(*this, N));
377  }
378 
379  /* Third slot. */
380  // TODO Allow multiple dimensions here
381  template <char i, char j, int Dim>
382  typename std::enable_if<
383  (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim),
387  T, Dim, i, j>>::type
388  operator()(const Index<i, Dim>, const Index<j, Dim>, const int N) const
389  {
390  using TensorExpr = Tensor3_numeral_2<
393  TensorExpr(*this, N));
394  }
395 
396  /* This is for expressions where a numeral is used for two slots, and
397  an Index for the other, yielding a Tensor1_Expr. */
398 
399  /* Index in first slot. */
400 
401  template <char i, int Dim>
402  typename std::enable_if<
403  (Tensor_Dim0 >= Dim),
406  T, Dim, i>>::type
407  operator()(const Index<i, Dim> index, const int N1, const int N2) const
408  {
409  using TensorExpr = Tensor3_numeral_12<
411  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
412  }
413 
414  /* Index in second slot. I use the same structures as for the Index
415  in the first slot since the tensor is symmetric on the first two
416  indices. */
417 
418  template <char i, int Dim>
419  typename std::enable_if<
420  (Tensor_Dim1 >= Dim),
423  T, Dim, i>>::type
424  operator()(const int N1, const Index<i, Dim> index, const int N2) const
425  {
426  using TensorExpr = Tensor3_numeral_12<
428  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
429  }
430 
431  /* Index in third slot. */
432 
433  template <char i, int Dim>
434  typename std::enable_if<
435  (Tensor_Dim2 >= Dim),
438  T, Dim, i>>::type
439  operator()(const int N1, const int N2, const Index<i, Dim> index) const
440  {
441  using TensorExpr = Tensor3_numeral_01<
443  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
444  }
445  };
446 }
447 
448 /// JSON compatible output for numbers. If T is something funky
449 /// (strings, arrays, etc.) then you are going to have a bad time.
450 
451 namespace FTensor
452 {
453  template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
454  std::ostream &Tensor3_ostream_row(
455  std::ostream &os,
457  const int &i, const int &j)
458  {
459  os << '[';
460  for(int k = 0; k + 1 < Tensor_Dim2; ++k)
461  {
462  os << t(i, j, k) << ',';
463  }
464  if(Tensor_Dim2 > 0)
465  {
466  os << t(i, j, Tensor_Dim2 - 1);
467  }
468  os << ']';
469  return os;
470  }
471 
472  template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
473  std::ostream &Tensor3_ostream_block(
474  std::ostream &os,
476  const int &i)
477  {
478  os << '[';
479  for(int j = 0; j + 1 < Tensor_Dim1; ++j)
480  {
481  FTensor::Tensor3_ostream_row(os, t, i, j);
482  os << ',';
483  }
484  if(Tensor_Dim1 > 0)
485  {
486  FTensor::Tensor3_ostream_row(os, t, i, Tensor_Dim1 - 1);
487  }
488  os << ']';
489  return os;
490  }
491 }
492 
493 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
494 std::ostream &
495 operator<<(std::ostream &os,
497 {
498  os << '[';
499  for(int i = 0; i + 1 < Tensor_Dim0; ++i)
500  {
502  os << ',';
503  }
504  if(Tensor_Dim0 > 0)
505  {
506  FTensor::Tensor3_ostream_block(os, t, Tensor_Dim0 - 1);
507  }
508  os << ']';
509  return os;
510 }
511 
512 /// JSON compatible input. It does not validate that separators are
513 /// actually braces '[' or commas ','. It also ignores errors from
514 /// missing trailing characters. So you could do something like
515 ///
516 /// Tensor3<double,2,2,2> t3_1;
517 /// std::stringstream ss(":::3:4:::7:8:::::11:12:::13:14");
518 /// ss >> t3_1;
519 
520 namespace FTensor
521 {
522  template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
523  std::istream &Tensor3_istream_row(
524  std::istream &is,
526  const int &i, const int &j)
527  {
528  char c;
529  is >> c;
530  for(int k = 0; k + 1 < Tensor_Dim2; ++k)
531  {
532  is >> t(i, j, k) >> c;
533  }
534  if(Tensor_Dim2 > 0)
535  {
536  is >> t(i, j, Tensor_Dim2 - 1);
537  }
538  is >> c;
539  return is;
540  }
541 
542  template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
543  std::istream &Tensor3_istream_block(
544  std::istream &is,
546  const int &i)
547  {
548  char c;
549  is >> c;
550  for(int j = 0; j + 1 < Tensor_Dim1; ++j)
551  {
552  Tensor3_istream_row(is, t, i, j);
553  is >> c;
554  }
555  if(Tensor_Dim1 > 0)
556  {
557  Tensor3_istream_row(is, t, i, Tensor_Dim1 - 1);
558  }
559  is >> c;
560  return is;
561  }
562 }
563 
564 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
565 std::istream &
566 operator>>(std::istream &is,
568 {
569  char c;
570  is >> c;
571  for(int i = 0; i + 1 < Tensor_Dim0; ++i)
572  {
574  is >> c;
575  }
576  if(Tensor_Dim0 > 0)
577  {
578  FTensor::Tensor3_istream_block(is, t, Tensor_Dim0 - 1);
579  }
580  is >> c;
581  return is;
582 }
T & operator()(const int N1, const int N2, const int N3)
T operator()(const int N1, const int N2, const int N3) const
Fully Antisymmetric Levi-Civita Tensor.
std::istream & Tensor3_istream_row(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i, const int &j)
std::istream & Tensor3_istream_block(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i)
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 & Tensor3_ostream_row(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i, const int &j)
std::istream & operator>>(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t)
std::ostream & operator<<(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t)
std::ostream & Tensor3_ostream_block(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i)
T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2]
const int N
Definition: speed_test.cpp:3