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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 //
9 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 //////////////////////////////////////////////////////////////////////////////////////
11
12
13 #ifndef ANTI_SYM_TENZOR_H
14 #define ANTI_SYM_TENZOR_H
15
16 /***************************************************************************
17 *
18 * The POOMA Framework
19 *
20 * This program was prepared by the Regents of the University of
21 * California at Los Alamos National Laboratory (the University) under
22 * Contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
23 * The University has certain rights in the program pursuant to the
24 * contract and the program should not be copied or distributed outside
25 * your organization. All rights in the program are reserved by the DOE
26 * and the University. Neither the U.S. Government nor the University
27 * makes any warranty, express or implied, or assumes any liability or
28 * responsibility for the use of this software
29 *
30 * Visit http://www.acl.lanl.gov/POOMA for more details
31 *
32 ***************************************************************************/
33
34
35 // include files
36 #include "PETE/PETE.h"
37 #include "OhmmsPETE/OhmmsTinyMeta.h"
38 #include "OhmmsPETE/Tensor.h"
39
40 #include <iostream>
41
42 namespace qmcplusplus
43 {
44 //////////////////////////////////////////////////////////////////////
45 //
46 // Definition of class AntiSymTensor.
47 //
48 //////////////////////////////////////////////////////////////////////
49
50 //
51 // | O -x10 -x20 |
52 // | x10 0 -x21 |
53 // | x20 x21 0 |
54 //
55
56 template<class T, unsigned D>
57 class AntiSymTensor
58 {
59 public:
60 typedef T Type_t;
61 enum
62 {
63 ElemDim = 2
64 };
65 enum
66 {
67 Size = D * (D - 1) / 2
68 };
69
70 // Default Constructor
AntiSymTensor()71 AntiSymTensor() { OTAssign<AntiSymTensor<T, D>, T, OpAssign>::apply(*this, T(0), OpAssign()); }
72
73 // A noninitializing ctor.
74 class DontInitialize
75 {};
AntiSymTensor(DontInitialize)76 AntiSymTensor(DontInitialize) {}
77
78 // Construct an AntiSymTensor from a single T.
79 // This doubles as the 2D AntiSymTensor initializer.
AntiSymTensor(const T & x00)80 AntiSymTensor(const T& x00) { OTAssign<AntiSymTensor<T, D>, T, OpAssign>::apply(*this, x00, OpAssign()); }
81 // construct a 3D AntiSymTensor
AntiSymTensor(const T & x10,const T & x20,const T & x21)82 AntiSymTensor(const T& x10, const T& x20, const T& x21)
83 {
84 X[0] = x10;
85 X[1] = x20;
86 X[2] = x21;
87 }
88
89 // Copy Constructor
AntiSymTensor(const AntiSymTensor<T,D> & rhs)90 AntiSymTensor(const AntiSymTensor<T, D>& rhs)
91 {
92 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T, D>, OpAssign>::apply(*this, rhs, OpAssign());
93 }
94
95 // Construct from a Tensor.
96 // Extract the antisymmetric part.
AntiSymTensor(const Tensor<T,D> & t)97 AntiSymTensor(const Tensor<T, D>& t)
98 {
99 for (int i = 1; i < D; ++i)
100 {
101 for (int j = 0; j < i; ++j)
102 (*this)[((i - 1) * i / 2) + j] = (t(i, j) - t(j, i)) * 0.5;
103 }
104 }
105
~AntiSymTensor()106 ~AntiSymTensor(){};
107
108 // assignment operators
109 const AntiSymTensor<T, D>& operator=(const AntiSymTensor<T, D>& rhs)
110 {
111 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T, D>, OpAssign>::apply(*this, rhs, OpAssign());
112 return *this;
113 }
114 template<class T1>
115 const AntiSymTensor<T, D>& operator=(const AntiSymTensor<T1, D>& rhs)
116 {
117 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T1, D>, OpAssign>::apply(*this, rhs, OpAssign());
118 return *this;
119 }
120 const AntiSymTensor<T, D>& operator=(const T& rhs)
121 {
122 OTAssign<AntiSymTensor<T, D>, T, OpAssign>::apply(*this, rhs, OpAssign());
123 return *this;
124 }
125
126 // accumulation operators
127 template<class T1>
128 AntiSymTensor<T, D>& operator+=(const AntiSymTensor<T1, D>& rhs)
129 {
130 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T1, D>, OpAddAssign>::apply(*this, rhs, OpAddAssign());
131 return *this;
132 }
133
134 template<class T1>
135 AntiSymTensor<T, D>& operator-=(const AntiSymTensor<T1, D>& rhs)
136 {
137 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T1, D>, OpSubtractAssign>::apply(*this, rhs, OpSubtractAssign());
138 return *this;
139 }
140
141 template<class T1>
142 AntiSymTensor<T, D>& operator*=(const AntiSymTensor<T1, D>& rhs)
143 {
144 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T1, D>, OpMultiplyAssign>::apply(*this, rhs, OpMultiplyAssign());
145 return *this;
146 }
147 AntiSymTensor<T, D>& operator*=(const T& rhs)
148 {
149 OTAssign<AntiSymTensor<T, D>, T, OpMultiplyAssign>::apply(*this, rhs, OpMultiplyAssign());
150 return *this;
151 }
152
153 template<class T1>
154 AntiSymTensor<T, D>& operator/=(const AntiSymTensor<T1, D>& rhs)
155 {
156 OTAssign<AntiSymTensor<T, D>, AntiSymTensor<T1, D>, OpDivideAssign>::apply(*this, rhs, OpDivideAssign());
157 return *this;
158 }
159 AntiSymTensor<T, D>& operator/=(const T& rhs)
160 {
161 OTAssign<AntiSymTensor<T, D>, T, OpDivideAssign>::apply(*this, rhs, OpDivideAssign());
162 return *this;
163 }
164
165 // Methods
166
len(void)167 int len(void) const { return Size; }
size(void)168 int size(void) const { return sizeof(*this); }
get_Size(void)169 int get_Size(void) const { return Size; }
170
171 class AssignProxy
172 {
173 public:
AssignProxy(Type_t & elem,int where)174 AssignProxy(Type_t& elem, int where) : elem_m(elem), where_m(where) {}
AssignProxy(const AssignProxy & model)175 AssignProxy(const AssignProxy& model) : elem_m(model.elem_m), where_m(where_m) {}
176 const AssignProxy& operator=(const AssignProxy& a)
177 {
178 PAssert(where_m != 0 || a.elem_m == -a.elem_m);
179 elem_m = where_m < 0 ? -a.elem_m : a.elem_m;
180 return *this;
181 }
182 const AssignProxy& operator=(const Type_t& e)
183 {
184 PAssert(where_m != 0 || e == -e);
185 elem_m = where_m < 0 ? -e : e;
186 return *this;
187 }
188
Type_t()189 operator Type_t() const { return (where_m < 0 ? -elem_m : elem_m); }
190
191 private:
192 //mutable Type_t &elem_m;
193 //mutable int where_m;
194 Type_t& elem_m;
195 int where_m;
196 };
197
198 // Operators
199
operator()200 Type_t operator()(unsigned int i, unsigned int j) const
201 {
202 if (i == j)
203 return T(0.0);
204 else if (i < j)
205 return -X[((j - 1) * j / 2) + i];
206 else
207 return X[((i - 1) * i / 2) + j];
208 }
209
210 // Type_t operator()( std::pair<int,int> a) const {
211 // int i = a.first;
212 // int j = a.second;
213 // return (*this)(i, j);
214 // }
215
operator()216 AssignProxy operator()(unsigned int i, unsigned int j)
217 {
218 if (i == j)
219 return AssignProxy(AntiSymTensor<T, D>::Zero, 0);
220 else
221 {
222 int lo = i < j ? i : j;
223 int hi = i > j ? i : j;
224 return AssignProxy(X[((hi - 1) * hi / 2) + lo], i - j);
225 }
226 }
227
228 // AssignProxy operator()(std::pair<int,int> a) {
229 // int i = a.first;
230 // int j = a.second;
231 // return (*this)(i, j);
232 // }
233
234 Type_t& operator[](unsigned int i)
235 {
236 PAssert(i < Size);
237 return X[i];
238 }
239
240 Type_t operator[](unsigned int i) const
241 {
242 PAssert(i < Size);
243 return X[i];
244 }
245
246 // These are the same as operator[] but with () instead:
247
operator()248 Type_t& operator()(unsigned int i)
249 {
250 PAssert(i < Size);
251 return X[i];
252 }
253
operator()254 Type_t operator()(unsigned int i) const
255 {
256 PAssert(i < Size);
257 return X[i];
258 }
259
260 // //----------------------------------------------------------------------
261 // // Comparison operators.
262 // bool operator==(const AntiSymTensor<T,D>& that) const {
263 // return TSV_MetaCompareArrays<T,T,D*(D-1)/2>::apply(X,that.X);
264 // }
265 // bool operator!=(const AntiSymTensor<T,D>& that) const {
266 // return !(*this == that);
267 // }
268
269 // //----------------------------------------------------------------------
270 // // parallel communication
271 // Message& putMessage(Message& m) const {
272 // m.setCopy(true);
273 // ::putMessage(m, X, X + Size);
274 // return m;
275 // }
276
277 // Message& getMessage(Message& m) {
278 // ::getMessage(m, X, X + Size);
279 // return m;
280 // }
281
282 private:
283 friend class AssignProxy;
284
285 // The elements themselves.
286 T X[Size];
287
288 // A place to store a zero element.
289 // We need to return a reference to this
290 // for the diagonal element.
291 static T Zero;
292 };
293
294
295 ///////////////////////////////////////////////////////////////////////////
296 // Specialization for 1D -- this is basically just the zero tensor
297 ///////////////////////////////////////////////////////////////////////////
298
299 template<class T>
300 class AntiSymTensor<T, 1>
301 {
302 public:
303 typedef T Type_t;
304 enum
305 {
306 ElemDim = 2
307 };
308
309 // Default Constructor
AntiSymTensor()310 AntiSymTensor() {}
311
312 // Copy Constructor
AntiSymTensor(const AntiSymTensor<T,1> &)313 AntiSymTensor(const AntiSymTensor<T, 1>&) {}
314
~AntiSymTensor()315 ~AntiSymTensor() {}
316
317 // assignment operators
318 const AntiSymTensor<T, 1>& operator=(const AntiSymTensor<T, 1>&) { return *this; }
319 template<class T1>
320 const AntiSymTensor<T, 1>& operator=(const AntiSymTensor<T1, 1>&)
321 {
322 return *this;
323 }
324 const AntiSymTensor<T, 1>& operator=(const T& rhs) { return *this; }
325
326 // accumulation operators
327 template<class T1>
328 AntiSymTensor<T, 1>& operator+=(const AntiSymTensor<T1, 1>&)
329 {
330 return *this;
331 }
332
333 template<class T1>
334 AntiSymTensor<T, 1>& operator-=(const AntiSymTensor<T1, 1>&)
335 {
336 return *this;
337 }
338
339 template<class T1>
340 AntiSymTensor<T, 1>& operator*=(const AntiSymTensor<T1, 1>&)
341 {
342 return *this;
343 }
344 AntiSymTensor<T, 1>& operator*=(const T&) { return *this; }
345
346 template<class T1>
347 AntiSymTensor<T, 1>& operator/=(const AntiSymTensor<T1, 1>&)
348 {
349 return *this;
350 }
351 AntiSymTensor<T, 1>& operator/=(const T&) { return *this; }
352
353 // Methods
354
len(void)355 int len(void) const { return Size; }
size(void)356 int size(void) const { return sizeof(*this); }
get_Size(void)357 int get_Size(void) const { return Size; }
358
359 class AssignProxy
360 {
361 public:
AssignProxy(Type_t & elem,int where)362 AssignProxy(Type_t& elem, int where) : elem_m(elem), where_m(where) {}
AssignProxy(const AssignProxy & model)363 AssignProxy(const AssignProxy& model) : elem_m(model.elem_m), where_m(where_m) {}
364 const AssignProxy& operator=(const AssignProxy& a)
365 {
366 PAssert(where_m != 0 || a.elem_m == -a.elem_m);
367 elem_m = where_m < 0 ? -a.elem_m : a.elem_m;
368 return *this;
369 }
370 const AssignProxy& operator=(const Type_t& e)
371 {
372 PAssert(where_m != 0 || e == -e);
373 elem_m = where_m < 0 ? -e : e;
374 return *this;
375 }
376
Type_t()377 operator Type_t() const { return (where_m < 0 ? -elem_m : elem_m); }
378
379 private:
380 //mutable Type_t &elem_m;
381 //mutable int where_m;
382 Type_t& elem_m;
383 int where_m;
384 };
385
386 // Operators
387
operator()388 Type_t operator()(unsigned int i, unsigned int j) const
389 {
390 PAssert(i == j);
391 return T(0.0);
392 }
393
394 // Type_t operator()( std::pair<int,int> a) const {
395 // int i = a.first;
396 // int j = a.second;
397 // return (*this)(i, j);
398 // }
399
operator()400 AssignProxy operator()(unsigned int i, unsigned int j)
401 {
402 PAssert(i == j);
403 return AssignProxy(AntiSymTensor<T, 1>::Zero, 0);
404 }
405
406 // AssignProxy operator()(std::pair<int,int> a) {
407 // int i = a.first;
408 // int j = a.second;
409 // return (*this)(i, j);
410 // }
411
412 Type_t operator[](unsigned int i) const
413 {
414 PAssert(i == 0);
415 return AntiSymTensor<T, 1>::Zero;
416 }
417
418 // These are the same as operator[] but with () instead:
419
operator()420 Type_t operator()(unsigned int i) const
421 {
422 PAssert(i == 0);
423 return AntiSymTensor<T, 1>::Zero;
424 }
425
426 //----------------------------------------------------------------------
427 // Comparison operators.
428 bool operator==(const AntiSymTensor<T, 1>& that) const { return true; }
429 bool operator!=(const AntiSymTensor<T, 1>& that) const { return !(*this == that); }
430
431 //----------------------------------------------------------------------
432 // parallel communication
433 // Message& putMessage(Message& m) const {
434 // m.setCopy(true);
435 // m.put(AntiSymTensor<T,1>::Zero);
436 // return m;
437 // }
438
439 // Message& getMessage(Message& m) {
440 // T zero;
441 // m.get(zero);
442 // return m;
443 // }
444
445 private:
446 friend class AssignProxy;
447
448 // The number of elements.
449 enum
450 {
451 Size = 0
452 };
453
454 // A place to store a zero element.
455 // We need to return a reference to this
456 // for the diagonal element.
457 static T Zero;
458 };
459
460
461 //////////////////////////////////////////////////////////////////////
462 //
463 // Free functions
464 //
465 //////////////////////////////////////////////////////////////////////
466
467 template<class T, unsigned D>
trace(const AntiSymTensor<T,D> &)468 T trace(const AntiSymTensor<T, D>&)
469 {
470 return T(0.0);
471 }
472
473 template<class T, unsigned D>
transpose(const AntiSymTensor<T,D> & rhs)474 AntiSymTensor<T, D> transpose(const AntiSymTensor<T, D>& rhs)
475 {
476 return -rhs;
477 }
478 //////////////////////////////////////////////////////////////////////
479 //
480 // Unary Operators
481 //
482 //////////////////////////////////////////////////////////////////////
483
484 //----------------------------------------------------------------------
485 // unary operator-
486 // template<class T, unsigned D>
487 // inline AntiSymTensor<T,D> operator-(const AntiSymTensor<T,D> &op)
488 // {
489 // return MetaUnary< AntiSymTensor<T,D> , OpUnaryMinus > :: apply(op,OpUnaryMinus());
490 // }
491
492 //----------------------------------------------------------------------
493 // unary operator+
494 template<class T, unsigned D>
495 inline const AntiSymTensor<T, D>& operator+(const AntiSymTensor<T, D>& op)
496 {
497 return op;
498 }
499
500 //////////////////////////////////////////////////////////////////////
501 //
502 // Binary Operators
503 //
504 //////////////////////////////////////////////////////////////////////
505
506 //
507 // Elementwise operators.
508 //
509
510 OHMMS_META_BINARY_OPERATORS(AntiSymTensor, operator+, OpAdd)
511 OHMMS_META_BINARY_OPERATORS(AntiSymTensor, operator-, OpSubtract)
OHMMS_META_BINARY_OPERATORS(AntiSymTensor,operator *,OpMultiply)512 OHMMS_META_BINARY_OPERATORS(AntiSymTensor, operator*, OpMultiply)
513 OHMMS_META_BINARY_OPERATORS(AntiSymTensor, operator/, OpDivide)
514
515 //----------------------------------------------------------------------
516 // dot products.
517 //----------------------------------------------------------------------
518
519 template<class T1, class T2, unsigned D>
520 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const AntiSymTensor<T1, D>& lhs,
521 const AntiSymTensor<T2, D>& rhs)
522 {
523 return OTDot<AntiSymTensor<T1, D>, AntiSymTensor<T2, D>>::apply(lhs, rhs);
524 }
525
526 template<class T1, class T2, unsigned D>
dot(const AntiSymTensor<T1,D> & lhs,const Tensor<T2,D> & rhs)527 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const AntiSymTensor<T1, D>& lhs,
528 const Tensor<T2, D>& rhs)
529 {
530 return OTDot<Tensor<T1, D>, Tensor<T2, D>>::apply(Tensor<T1, D>(lhs), rhs);
531 }
532
533 template<class T1, class T2, unsigned D>
dot(const Tensor<T1,D> & lhs,const AntiSymTensor<T2,D> & rhs)534 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const Tensor<T1, D>& lhs,
535 const AntiSymTensor<T2, D>& rhs)
536 {
537 return OTDot<Tensor<T1, D>, Tensor<T2, D>>::apply(lhs, Tensor<T2, D>(rhs));
538 }
539
540 template<class T1, class T2, unsigned D>
dot(const AntiSymTensor<T1,D> & lhs,const SymTensor<T2,D> & rhs)541 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const AntiSymTensor<T1, D>& lhs,
542 const SymTensor<T2, D>& rhs)
543 {
544 return OTDot<Tensor<T1, D>, Tensor<T2, D>>::apply(Tensor<T1, D>(lhs), Tensor<T2, D>(rhs));
545 }
546
547 template<class T1, class T2, unsigned D>
dot(const SymTensor<T1,D> & lhs,const AntiSymTensor<T2,D> & rhs)548 inline Tensor<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const SymTensor<T1, D>& lhs,
549 const AntiSymTensor<T2, D>& rhs)
550 {
551 return OTDot<Tensor<T1, D>, Tensor<T2, D>>::apply(Tensor<T1, D>(lhs), Tensor<T2, D>(rhs));
552 }
553
554 template<class T1, class T2, unsigned D>
dot(const TinyVector<T1,D> & lhs,const AntiSymTensor<T2,D> & rhs)555 inline TinyVector<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const TinyVector<T1, D>& lhs,
556 const AntiSymTensor<T2, D>& rhs)
557 {
558 return OTDot<TinyVector<T1, D>, AntiSymTensor<T2, D>>::apply(lhs, rhs);
559 }
560
561 template<class T1, class T2, unsigned D>
dot(const AntiSymTensor<T1,D> & lhs,const TinyVector<T2,D> & rhs)562 inline TinyVector<typename BinaryReturn<T1, T2, OpMultiply>::Type_t, D> dot(const AntiSymTensor<T1, D>& lhs,
563 const TinyVector<T2, D>& rhs)
564 {
565 return OTDot<AntiSymTensor<T1, D>, TinyVector<T2, D>>::apply(lhs, rhs);
566 }
567
568 //----------------------------------------------------------------------
569 // double dot products.
570 //----------------------------------------------------------------------
571
572 //template < class T1, class T2, unsigned D >
573 //inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
574 //dotdot(const AntiSymTensor<T1,D> &lhs, const AntiSymTensor<T2,D> &rhs)
575 //{
576 // return MetaDotDot< AntiSymTensor<T1,D> , AntiSymTensor<T2,D> > ::
577 // apply(lhs,rhs);
578 //}
579
580 //template < class T1, class T2, unsigned D >
581 //inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
582 //dotdot(const AntiSymTensor<T1,D> &lhs, const Tensor<T2,D> &rhs)
583 //{
584 // return MetaDotDot< Tensor<T1,D> , Tensor<T2,D> > ::
585 // apply(Tensor<T1,D>(lhs),rhs);
586 //}
587
588 //template < class T1, class T2, unsigned D >
589 //inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
590 //dotdot(const Tensor<T1,D> &lhs, const AntiSymTensor<T2,D> &rhs)
591 //{
592 // return MetaDotDot< Tensor<T1,D> , Tensor<T2,D> > ::
593 // apply(lhs,Tensor<T2,D>(rhs));
594 //}
595
596 //template < class T1, class T2, unsigned D >
597 //inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
598 //dotdot(const AntiSymTensor<T1,D> &lhs, const SymTensor<T2,D> &rhs)
599 //{
600 // return MetaDotDot< Tensor<T1,D> , Tensor<T2,D> > ::
601 // apply(Tensor<T1,D>(lhs),Tensor<T2,D>(rhs));
602 //}
603
604 //template < class T1, class T2, unsigned D >
605 //inline typename BinaryReturn<T1,T2,OpMultiply>::Type_t
606 //dotdot(const SymTensor<T1,D> &lhs, const AntiSymTensor<T2,D> &rhs)
607 //{
608 // return MetaDotDot< Tensor<T1,D> , Tensor<T2,D> > ::
609 // apply(Tensor<T1,D>(lhs),Tensor<T2,D>(rhs));
610 //}
611
612 //----------------------------------------------------------------------
613 // I/O
614 template<class T, unsigned D>
615 std::ostream& operator<<(std::ostream& out, const AntiSymTensor<T, D>& rhs)
616 {
617 if (D >= 1)
618 {
619 for (int i = 0; i < D; i++)
620 {
621 for (int j = 0; j < D - 1; j++)
622 {
623 out << rhs(i, j) << " ";
624 }
625 out << rhs(i, D - 1) << " ";
626 if (i < D - 1)
627 out << std::endl;
628 }
629 }
630 else
631 {
632 out << " " << rhs(0, 0) << " ";
633 }
634 return out;
635 }
636
637 //////////////////////////////////////////////////////////////////////
638
639 template<class T, unsigned int D>
640 T AntiSymTensor<T, D>::Zero = 0;
641
642 } // namespace qmcplusplus
643 #endif // ANTI_SYM_TENZOR_H
644