1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
10 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 //
12 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
13 //////////////////////////////////////////////////////////////////////////////////////
14
15
16 #ifndef OHMMS_TENSOR_H
17 #define OHMMS_TENSOR_H
18
19 #include "PETE/PETE.h"
20 #include "OhmmsPETE/OhmmsTinyMeta.h"
21 /***************************************************************************
22 *
23 * The POOMA Framework
24 *
25 * This program was prepared by the Regents of the University of
26 * California at Los Alamos National Laboratory (the University) under
27 * Contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
28 * The University has certain rights in the program pursuant to the
29 * contract and the program should not be copied or distributed outside
30 * your organization. All rights in the program are reserved by the DOE
31 * and the University. Neither the U.S. Government nor the University
32 * makes any warranty, express or implied, or assumes any liability or
33 * responsibility for the use of this software
34 *
35 * Visit http://www.acl.lanl.gov/POOMA for more details
36 *
37 ***************************************************************************/
38
39
40 namespace qmcplusplus
41 {
42 // forward declarations
43 template<class T, unsigned D>
44 class SymTensor;
45 template<class T, unsigned D>
46 class AntiSymTensor;
47
48
49 /** Tensor<T,D> class for D by D tensor
50 *
51 * @tparam T datatype
52 * @tparm D dimension
53 */
54 template<class T, unsigned D>
55 class Tensor
56 {
57 public:
58 typedef T Type_t;
59 enum
60 {
61 ElemDim = 2
62 };
63 enum
64 {
65 Size = D * D
66 };
67
68 // Default Constructor
Tensor()69 inline Tensor()
70 {
71 for (size_t d = 0; d < Size; ++d)
72 X[d] = T(0);
73 }
74
75 // A noninitializing ctor.
76 class DontInitialize
77 {};
Tensor(DontInitialize)78 inline Tensor(DontInitialize) {}
79
80 // Copy Constructor
81 inline Tensor(const Tensor& rhs) = default;
82
83 template<class T1>
Tensor(const Tensor<T1,D> & rhs)84 inline Tensor(const Tensor<T1, D>& rhs)
85 {
86 for (size_t d = 0; d < Size; ++d)
87 X[d] = rhs[d];
88 }
89
90 // constructor from a single T
Tensor(const T & x00)91 inline Tensor(const T& x00)
92 {
93 for (size_t d = 0; d < Size; ++d)
94 X[d] = x00;
95 }
96
97 // constructors for fixed dimension
Tensor(const T & x00,const T & x10,const T & x01,const T & x11)98 Tensor(const T& x00, const T& x10, const T& x01, const T& x11)
99 {
100 X[0] = x00;
101 X[1] = x10;
102 X[2] = x01;
103 X[3] = x11;
104 }
Tensor(const T & x00,const T & x10,const T & x20,const T & x01,const T & x11,const T & x21,const T & x02,const T & x12,const T & x22)105 Tensor(const T& x00,
106 const T& x10,
107 const T& x20,
108 const T& x01,
109 const T& x11,
110 const T& x21,
111 const T& x02,
112 const T& x12,
113 const T& x22)
114 {
115 X[0] = x00;
116 X[1] = x10;
117 X[2] = x20;
118 X[3] = x01;
119 X[4] = x11;
120 X[5] = x21;
121 X[6] = x02;
122 X[7] = x12;
123 X[8] = x22;
124 }
125
126 //constructor from SymTensor
127 Tensor(const SymTensor<T, D>&);
128
129 // constructor from AntiSymTensor
130 Tensor(const AntiSymTensor<T, D>&);
131
132 // destructor
~Tensor()133 ~Tensor(){};
134
135 // assignment operators
136 inline Tensor& operator=(const Tensor& rhs) = default;
137
138 template<class T1>
139 inline Tensor<T, D>& operator=(const Tensor<T1, D>& rhs)
140 {
141 OTAssign<Tensor<T, D>, Tensor<T1, D>, OpAssign>::apply(*this, rhs, OpAssign());
142 return *this;
143 }
144 inline Tensor<T, D>& operator=(const T& rhs)
145 {
146 OTAssign<Tensor<T, D>, T, OpAssign>::apply(*this, rhs, OpAssign());
147 return *this;
148 }
149
150 // accumulation operators
151 template<class T1>
152 inline Tensor<T, D>& operator+=(const Tensor<T1, D>& rhs)
153 {
154 OTAssign<Tensor<T, D>, Tensor<T1, D>, OpAddAssign>::apply(*this, rhs, OpAddAssign());
155 return *this;
156 }
157 inline Tensor<T, D>& operator+=(const T& rhs)
158 {
159 OTAssign<Tensor<T, D>, T, OpAddAssign>::apply(*this, rhs, OpAddAssign());
160 return *this;
161 }
162
163 template<class T1>
164 inline Tensor<T, D>& operator-=(const Tensor<T1, D>& rhs)
165 {
166 OTAssign<Tensor<T, D>, Tensor<T1, D>, OpSubtractAssign>::apply(*this, rhs, OpSubtractAssign());
167 return *this;
168 }
169
170 inline Tensor<T, D>& operator-=(const T& rhs)
171 {
172 OTAssign<Tensor<T, D>, T, OpSubtractAssign>::apply(*this, rhs, OpSubtractAssign());
173 return *this;
174 }
175
176 template<class T1>
177 inline Tensor<T, D>& operator*=(const Tensor<T1, D>& rhs)
178 {
179 OTAssign<Tensor<T, D>, Tensor<T1, D>, OpMultiplyAssign>::apply(*this, rhs, OpMultiplyAssign());
180 return *this;
181 }
182
183 inline Tensor<T, D>& operator*=(const T& rhs)
184 {
185 OTAssign<Tensor<T, D>, T, OpMultiplyAssign>::apply(*this, rhs, OpMultiplyAssign());
186 return *this;
187 }
188
189 template<class T1>
190 inline Tensor<T, D>& operator/=(const Tensor<T1, D>& rhs)
191 {
192 OTAssign<Tensor<T, D>, Tensor<T1, D>, OpDivideAssign>::apply(*this, rhs, OpDivideAssign());
193 return *this;
194 }
195
196 inline Tensor<T, D>& operator/=(const T& rhs)
197 {
198 OTAssign<Tensor<T, D>, T, OpDivideAssign>::apply(*this, rhs, OpDivideAssign());
199 return *this;
200 }
201
202 // Methods
203
diagonal(const T & rhs)204 inline void diagonal(const T& rhs)
205 {
206 for (int i = 0; i < D; i++)
207 (*this)(i, i) = rhs;
208 }
209
add2diagonal(T rhs)210 inline void add2diagonal(T rhs)
211 {
212 for (int i = 0; i < D; i++)
213 (*this)(i, i) += rhs;
214 }
215
216 ///return the size
len()217 inline int len() const { return Size; }
218 ///return the size
size()219 inline int size() const { return Size; }
220
221 /** return the i-th value or assign
222 * @param i index [0,D*D)
223 */
224 inline Type_t& operator[](unsigned int i) { return X[i]; }
225
226 /** return the i-th value
227 * @param i index [0,D*D)
228 */
229 inline Type_t operator[](unsigned int i) const { return X[i]; }
230
231 //TJW: add these 12/16/97 to help with NegReflectAndZeroFace BC:
232 // These are the same as operator[] but with () instead:
operator()233 inline Type_t& operator()(unsigned int i) { return X[i]; }
234
operator()235 inline Type_t operator()(unsigned int i) const { return X[i]; }
236 //TJW.
237
238 /** return the (i,j) component
239 * @param i index [0,D)
240 * @param j index [0,D)
241 */
operator()242 inline Type_t operator()(unsigned int i, unsigned int j) const { return X[i * D + j]; }
243
244 /** return/assign the (i,j) component
245 * @param i index [0,D)
246 * @param j index [0,D)
247 */
operator()248 inline Type_t& operator()(unsigned int i, unsigned int j) { return X[i * D + j]; }
249
getRow(unsigned int i)250 inline TinyVector<T, D> getRow(unsigned int i)
251 {
252 TinyVector<T, D> res;
253 for (int j = 0; j < D; j++)
254 res[j] = X[i * D + j];
255 return res;
256 }
257
getColumn(unsigned int i)258 inline TinyVector<T, D> getColumn(unsigned int i)
259 {
260 TinyVector<T, D> res;
261 for (int j = 0; j < D; j++)
262 res[j] = X[j * D + i];
263 return res;
264 }
265
data()266 inline Type_t* data() { return X; }
data()267 inline const Type_t* data() const { return X; }
begin()268 inline Type_t* begin() { return X; }
begin()269 inline const Type_t* begin() const { return X; }
end()270 inline Type_t* end() { return X + Size; }
end()271 inline const Type_t* end() const { return X + Size; }
272
273 // Removed operator using std::pair
274 // Type_t operator()(const std::pair<int,int> i) const {
275 // PAssert ( (i.first>=0) && (i.second>=0) && (i.first<D) && (i.second<D) );
276 // return (*this)(i.first,i.second);
277 // }
278
279 // Type_t& operator()(const std::pair<int,int> i) {
280 // PAssert ( (i.first>=0) && (i.second>=0) && (i.first<D) && (i.second<D) );
281 // return (*this)(i.first,i.second);
282 // }
283
284
285 //----------------------------------------------------------------------
286 // Comparison operators.
287 // bool operator==(const Tensor<T,D>& that) const {
288 // return MetaCompareArrays<T,T,D*D>::apply(X,that.X);
289 // }
290 // bool operator!=(const Tensor<T,D>& that) const {
291 // return !(*this == that);
292 // }
293
294 //----------------------------------------------------------------------
295 // parallel communication
296 // Message& putMessage(Message& m) const {
297 // m.setCopy(true);
298 // const T *p = X;
299 // ::putMessage(m, p, p + D*D);
300 // return m;
301 // }
302
303 // Message& getMessage(Message& m) {
304 // T *p = X;
305 // ::getMessage(m, p, p + D*D);
306 // return m;
307 // }
308
309 private:
310 // The elements themselves.
311 T X[Size];
312 };
313
314
315 //////////////////////////////////////////////////////////////////////
316 //
317 // Free functions
318 //
319 //////////////////////////////////////////////////////////////////////
320
321 /** trace \f$ result = \sum_k rhs(k,k)\f$
322 * @param rhs a tensor
323 */
324 template<class T, unsigned D>
trace(const Tensor<T,D> & rhs)325 inline T trace(const Tensor<T, D>& rhs)
326 {
327 T result = 0.0;
328 for (int i = 0; i < D; i++)
329 result += rhs(i, i);
330 return result;
331 }
332
333 /** transpose a tensor
334 */
335 template<class T, unsigned D>
transpose(const Tensor<T,D> & rhs)336 inline Tensor<T, D> transpose(const Tensor<T, D>& rhs)
337 {
338 Tensor<T, D> result; // = Tensor<T,D>::DontInitialize();
339 for (int j = 0; j < D; j++)
340 for (int i = 0; i < D; i++)
341 result(i, j) = rhs(j, i);
342 return result;
343 }
344
345 /** Tr(a*b), \f$ \sum_i\sum_j a(i,j)*b(j,i) \f$
346 */
347 template<class T1, class T2, unsigned D>
trace(const Tensor<T1,D> & a,const Tensor<T2,D> & b)348 inline T1 trace(const Tensor<T1, D>& a, const Tensor<T2, D>& b)
349 {
350 T1 result = 0.0;
351 for (int i = 0; i < D; i++)
352 for (int j = 0; j < D; j++)
353 result += a(i, j) * b(j, i);
354 return result;
355 }
356
357 /** Tr(a^t *b), \f$ \sum_i\sum_j a(i,j)*b(i,j) \f$
358 */
359 template<class T, unsigned D>
traceAtB(const Tensor<T,D> & a,const Tensor<T,D> & b)360 inline T traceAtB(const Tensor<T, D>& a, const Tensor<T, D>& b)
361 {
362 T result = 0.0;
363 for (int i = 0; i < D * D; i++)
364 result += a(i) * b(i);
365 return result;
366 }
367
368 /** Tr(a^t *b), \f$ \sum_i\sum_j a(i,j)*b(i,j) \f$
369 */
370 template<class T1, class T2, unsigned D>
traceAtB(const Tensor<T1,D> & a,const Tensor<T2,D> & b)371 inline typename BinaryReturn<T1, T2, OpMultiply>::Type_t traceAtB(const Tensor<T1, D>& a, const Tensor<T2, D>& b)
372 {
373 typedef typename BinaryReturn<T1, T2, OpMultiply>::Type_t T;
374 T result = 0.0;
375 for (int i = 0; i < D * D; i++)
376 result += a(i) * b(i);
377 return result;
378 }
379
380 //////////////////////////////////////////////////////////////////////
381 //
382 // Unary Operators
383 //
384 //////////////////////////////////////////////////////////////////////
385 //----------------------------------------------------------------------
386 // unary operator-
387 //template<class T, unsigned D>
388 //inline Tensor<T,D> operator-(const Tensor<T,D> &op)
389 //{
390 // return MetaUnary< Tensor<T,D> , OpUnaryMinus > :: apply(op,OpUnaryMinus());
391 //}
392 //----------------------------------------------------------------------
393 // unary operator+
394 //template<class T, unsigned D>
395 //inline const Tensor<T,D> &operator+(const Tensor<T,D> &op)
396 //{
397 // return op;
398 //}
399
400 /// Binary Operators
401 OHMMS_META_BINARY_OPERATORS(Tensor, operator+, OpAdd)
402 OHMMS_META_BINARY_OPERATORS(Tensor, operator-, OpSubtract)
OHMMS_META_BINARY_OPERATORS(Tensor,operator *,OpMultiply)403 OHMMS_META_BINARY_OPERATORS(Tensor, operator*, OpMultiply)
404 OHMMS_META_BINARY_OPERATORS(Tensor, operator/, OpDivide)
405 //TSV_ELEMENTWISE_OPERATOR(Tensor,Min,FnMin)
406 //TSV_ELEMENTWISE_OPERATOR(Tensor,Max,FnMax)
407
408 /** Tensor-Tensor dot product \f$result(i,j)=\sum_k lhs(i,k)*rhs(k,j)\f$
409 * @param lhs a tensor
410 * @param rhs a tensor
411 */
412 template<class T1, class T2, unsigned D>
413 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const Tensor<T1, D>& lhs,
414 const Tensor<T2, D>& rhs)
415 {
416 return OTDot<Tensor<T1, D>, Tensor<T2, D>>::apply(lhs, rhs);
417 }
418
419 /** Vector-Tensor dot product \f$result(i)=\sum_k lhs(k)*rhs(k,i)\f$
420 * @param lhs a vector
421 * @param rhs a tensor
422 */
423 template<class T1, class T2, unsigned D>
dot(const TinyVector<T1,D> & lhs,const Tensor<T2,D> & rhs)424 inline TinyVector<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const TinyVector<T1, D>& lhs,
425 const Tensor<T2, D>& rhs)
426 {
427 return OTDot<TinyVector<T1, D>, Tensor<T2, D>>::apply(lhs, rhs);
428 }
429
430 /** Tensor-Vector dot product \f$result(i)=\sum_k lhs(i,k)*rhs(k)\f$
431 * @param lhs a tensor
432 * @param rhs a vector
433 */
434 template<class T1, class T2, unsigned D>
dot(const Tensor<T1,D> & lhs,const TinyVector<T2,D> & rhs)435 inline TinyVector<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const Tensor<T1, D>& lhs,
436 const TinyVector<T2, D>& rhs)
437 {
438 return OTDot<Tensor<T1, D>, TinyVector<T2, D>>::apply(lhs, rhs);
439 }
440
441 //----------------------------------------------------------------------
442 // double dot product.
443 //----------------------------------------------------------------------
444
445 // template < class T1, class T2, unsigned D >
446 // inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
447 // dotdot(const Tensor<T1,D> &lhs, const Tensor<T2,D> &rhs)
448 // {
449 // return MetaDotDot< Tensor<T1,D> , Tensor<T2,D> > :: apply(lhs,rhs);
450 //}
451
452
453 //----------------------------------------------------------------------
454 // Outer product.
455 //----------------------------------------------------------------------
456
457 ///** Vector-vector outter product \f$ result(i,j)=v1(i)*v2(j)\f$
458 // * @param v1 a vector
459 // * @param v2 a vector
460 // */
461 // template<class T1, class T2, unsigned int D>
462 // inline Tensor<typename BinaryReturn<T1,T2,OpMultiply>::type,D >
463 // outerProduct(const TinyVector<T1,D>& v1, const TinyVector<T2,D>& v2)
464 // {
465 // typedef typename BinaryReturn<T1,T2,OpMultiply>::Type_t T0;
466 //// //#if (defined(POOMA_SGI_CC_721_TYPENAME_BUG) || defined(__MWERKS__))
467 //// //Tensor<T0,D> ret = Tensor<T0,D>::DontInitialize();
468 //// //#else
469 // Tensor<T0,D> ret = typename Tensor<T0,D>::DontInitialize();
470 //// //#endif // POOMA_SGI_CC_721_TYPENAME_BUG
471 // for (unsigned int i=0; i<D; ++i)
472 // for (unsigned int j=0; j<D; ++j)
473 // ret(i,j) = v1[i]*v2[j];
474 // return ret;
475 // }
476 //
477 // template<class T1, unsigned int D>
478 // inline Tensor<T1,D >
479 // outerProduct(const TinyVector<T1,D>& v1, const TinyVector<T1,D>& v2)
480 // {
481 // Tensor<T1,D> ret = typename Tensor<T1,D>::DontInitialize();
482 //// //#endif // POOMA_SGI_CC_721_TYPENAME_BUG
483 // for (unsigned int i=0; i<D; ++i)
484 // for (unsigned int j=0; j<D; ++j)
485 // ret(i,j) = v1[i]*v2[j];
486 // return ret;
487 // }
488 //
489
490
491 //----------------------------------------------------------------------
492 // I/O
493 template<class T, unsigned D>
494 std::ostream& operator<<(std::ostream& out, const Tensor<T, D>& rhs)
495 {
496 if (D >= 1)
497 {
498 for (int i = 0; i < D; i++)
499 {
500 for (int j = 0; j < D - 1; j++)
501 {
502 out << rhs(i, j) << " ";
503 }
504 out << rhs(i, D - 1) << " ";
505 if (i < D - 1)
506 out << std::endl;
507 }
508 }
509 else
510 {
511 out << " " << rhs(0, 0) << " ";
512 }
513 return out;
514 }
515
516 template<class T, unsigned D>
517 std::istream& operator>>(std::istream& is, Tensor<T, D>& rhs)
518 {
519 for (int i = 0; i < D * D; i++)
520 is >> rhs[i];
521 return is;
522 }
523
524 template<class T, unsigned D>
525 bool operator==(const Tensor<T, D>& lhs, const Tensor<T, D>& rhs)
526 {
527 for (int i = 0; i < D * D; i++)
528 if (lhs[i] != rhs[i])
529 return false;
530 return true;
531 }
532
533 template<class T, unsigned D>
534 bool operator!=(const Tensor<T, D>& lhs, const Tensor<T, D>& rhs)
535 {
536 return !(lhs == rhs);
537 }
538
539 } // namespace qmcplusplus
540 // include header files for SymTensor and AntiSymTensor in order
541 // to define constructors for Tensor using these types
542 #include "OhmmsPETE/SymTensor.h"
543 #include "OhmmsPETE/AntiSymTensor.h"
544
545 namespace qmcplusplus
546 {
547 template<class T, unsigned D>
Tensor(const SymTensor<T,D> & rhs)548 Tensor<T, D>::Tensor(const SymTensor<T, D>& rhs)
549 {
550 for (int i = 0; i < D; ++i)
551 for (int j = 0; j < D; ++j)
552 (*this)(i, j) = rhs(i, j);
553 }
554
555 template<class T, unsigned D>
Tensor(const AntiSymTensor<T,D> & rhs)556 Tensor<T, D>::Tensor(const AntiSymTensor<T, D>& rhs)
557 {
558 for (int i = 0; i < D; ++i)
559 for (int j = 0; j < D; ++j)
560 (*this)(i, j) = rhs(i, j);
561 }
562
563 } // namespace qmcplusplus
564
565 #endif // OHMMS_TENSOR_H
566