OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
IpplExpressions.h
Go to the documentation of this file.
1// -*- C++ -*-
2/***************************************************************************
3 *
4 * The IPPL Framework
5 *
6 *
7 * Visit http://people.web.psi.ch/adelmann/ for more details
8 *
9 ***************************************************************************/
10
11
13//
14// FILE NAME
15// IpplExpressions.h
16//
17// CREATED
18// July 11, 1997
19//
20// DESCRIPTION
21// This header file defines custom objects and operators necessary to use
22// expression templates in IPPL.
23//
25
26#ifndef IPPL_EXPRESSIONS_H
27#define IPPL_EXPRESSIONS_H
28
29
30// We need to construct a custom version of Reduction. We must define
31// this macro before including PETE.h.
32#define PETE_USER_REDUCTION_CODE \
33 R global_ret; \
34 reduce_masked(ret, global_ret, acc_op, 0 < n ); \
35 ret = global_ret;
36
37// include files
38#include "Message/Message.h"
40#include "PETE/PETE.h"
41
42
43//=========================================================================
44//
45// UNARY OPERATIONS
46//
47//=========================================================================
48
49// Abs is handled rather strangely. There appears to be two versions
50// that do the same thing: abs and Abs.
51
52PETE_DefineUnary(Abs, (0 < a ? a : -a), FnAbs)
53
54inline double
55PETE_apply(FnAbs, std::complex<double> a)
56{
57 return abs(a);
58}
59
60template<class T>
63{
65 (l.PETE_unwrap().MakeExpression());
66}
67
70PETE_DefineUnary(norm, (norm(a)), FnNorm)
78
79
80//=========================================================================
81//
82// BINARY OPERATIONS
83//
84//=========================================================================
85
86// define min/max for built-in scalar types
87#define PETE_DefineScalarMinMax(Sca) \
88inline Sca \
89Min(const Sca& a, const Sca& b) \
90{ \
91 return (a<b ? a : b); \
92} \
93inline Sca \
94Max(const Sca& a, const Sca& b) \
95{ \
96 return (b<a ? a : b); \
97}
98
104
107
112
119
120//tjw: make sure kosher to have "cross" down there (added in)
121#define PETE_DefineIPPLScalar(Sca) \
122PETE_DefineBinaryWithScalars(Min, FnMin, Sca) \
123PETE_DefineBinaryWithScalars(Max, FnMax, Sca) \
124PETE_DefineBinaryWithScalars(dot, FnDot, Sca) \
125PETE_DefineBinaryWithScalars(dotdot, FnDotDot, Sca) \
126PETE_DefineBinaryWithScalars(outerProduct, FnOuterProduct, Sca) \
127PETE_DefineBinaryWithScalars(cross, FnDot, Sca) \
128PETE_DefineBinaryWithScalars(lt, OpLT, Sca) \
129PETE_DefineBinaryWithScalars(le, OpLE, Sca) \
130PETE_DefineBinaryWithScalars(gt, OpGT, Sca) \
131PETE_DefineBinaryWithScalars(ge, OpGE, Sca) \
132PETE_DefineBinaryWithScalars(eq, OpEQ, Sca) \
133PETE_DefineBinaryWithScalars(ne, OpNE, Sca)
134
140
144
145#undef PETE_DefineIPPLScalar
146
147// Now we need to provide special cases for Vektors, SymTenzors, and Tenzors
148// so we can remove PETE_Expr base class.
149
150#define PETE_DefineBinaryWithVSTScalars(Fun,Op,Sca) \
151template<class T1, unsigned Dim, class T2> \
152inline PETE_TBTree<Op, PETE_Scalar< Sca<T1, Dim> >, \
153 typename T2::PETE_Expr_t> \
154Fun(const Sca<T1, Dim> &l, const PETE_Expr<T2>& r) \
155{ \
156 typedef PETE_TBTree<Op, PETE_Scalar< Sca<T1, Dim> >, \
157 typename T2::PETE_Expr_t> ret; \
158 return ret(PETE_Scalar< Sca<T1, Dim> >(l), \
159 r.PETE_unwrap().MakeExpression()); \
160} \
161template<class T1, class T2, unsigned Dim> \
162inline PETE_TBTree<Op, typename T1::PETE_Expr_t, \
163 PETE_Scalar< Sca<T2, Dim> > > \
164Fun(const PETE_Expr<T1>& l, const Sca<T2, Dim> &r) \
165{ \
166 typedef PETE_TBTree<Op, typename T1::PETE_Expr_t, \
167 PETE_Scalar< Sca<T2, Dim> > > ret; \
168 return ret(l.PETE_unwrap().MakeExpression(), \
169 PETE_Scalar< Sca<T2, Dim> >(r)); \
170}
171
172#define PETE_DefineTrinaryWithVSTScalars(Fun, Op, Sca) \
173template<class Cond_t, class True_t, class T, unsigned Dim> \
174inline PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, \
175 typename True_t::PETE_Expr_t, PETE_Scalar< Sca<T, Dim> > > \
176Fun(const PETE_Expr<Cond_t>& c, const PETE_Expr<True_t>& t, \
177 const Sca<T, Dim> &f) \
178{ \
179 typedef PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, \
180 typename True_t::PETE_Expr_t, PETE_Scalar< Sca<T, Dim> > > ret; \
181 return ret(c.PETE_unwrap().MakeExpression(), \
182 t.PETE_unwrap().MakeExpression(), PETE_Scalar< Sca<T, Dim> >(f)); \
183} \
184template<class Cond_t, class T, unsigned Dim, class False_t> \
185inline PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, \
186 PETE_Scalar< Sca<T, Dim> >, typename False_t::PETE_Expr_t > \
187Fun(const PETE_Expr<Cond_t>& c, const Sca<T, Dim> &t, \
188 const PETE_Expr<False_t>& f) \
189{ \
190 typedef PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, PETE_Scalar< Sca<T, Dim> >, \
191 typename False_t::PETE_Expr_t > ret; \
192 return ret(c.PETE_unwrap().MakeExpression(), \
193 PETE_Scalar< Sca<T, Dim> >(t), f.PETE_unwrap().MakeExpression()); \
194} \
195template<class Cond_t, class T, unsigned Dim> \
196inline PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, \
197 PETE_Scalar< Sca<T, Dim> >, PETE_Scalar< Sca<T, Dim> > > \
198Fun(const PETE_Expr<Cond_t>& c, const Sca<T, Dim> &t, const Sca<T, Dim> &f) \
199{ \
200 typedef PETE_TTTree<Op, typename Cond_t::PETE_Expr_t, \
201 PETE_Scalar< Sca<T, Dim> >, PETE_Scalar< Sca<T, Dim> > > ret; \
202 return ret(c.PETE_unwrap().MakeExpression(), \
203 PETE_Scalar< Sca<T, Dim> >(t), PETE_Scalar< Sca<T, Dim> >(f)); \
204}
205
206
207//tjw: make sure kosher to have "cross" down there (added in)
208#define PETE_DefineVSTScalar(Sca) \
209PETE_DefineBinaryWithVSTScalars(operator+, OpAdd, Sca) \
210PETE_DefineBinaryWithVSTScalars(operator-, OpSubtract, Sca) \
211PETE_DefineBinaryWithVSTScalars(operator*, OpMultipply, Sca) \
212PETE_DefineBinaryWithVSTScalars(operator/, OpDivide, Sca) \
213PETE_DefineBinaryWithVSTScalars(operator%, OpMod, Sca) \
214PETE_DefineBinaryWithVSTScalars(operator<, OpLT, Sca) \
215PETE_DefineBinaryWithVSTScalars(operator<=, OpLE, Sca) \
216PETE_DefineBinaryWithVSTScalars(operator>, OpGT, Sca) \
217PETE_DefineBinaryWithVSTScalars(operator>=, OpGE, Sca) \
218PETE_DefineBinaryWithVSTScalars(operator==, OpEQ, Sca) \
219PETE_DefineBinaryWithVSTScalars(operator!=, OpNE, Sca) \
220PETE_DefineBinaryWithVSTScalars(operator&&, OpAnd, Sca) \
221PETE_DefineBinaryWithVSTScalars(operator||, OpOr, Sca) \
222PETE_DefineBinaryWithVSTScalars(operator&, OpBitwiseAnd, Sca) \
223PETE_DefineBinaryWithVSTScalars(operator|, OpBitwiseOr, Sca) \
224PETE_DefineBinaryWithVSTScalars(operator^, OpBitwiseXor, Sca) \
225PETE_DefineBinaryWithVSTScalars(copysign, FnCopysign, Sca) \
226PETE_DefineBinaryWithVSTScalars(ldexp, FnLdexp, Sca) \
227PETE_DefineBinaryWithVSTScalars(pow, FnPow, Sca) \
228PETE_DefineBinaryWithVSTScalars(fmod, FnFmod, Sca) \
229PETE_DefineBinaryWithVSTScalars(atan2, FnArcTan2, Sca) \
230PETE_DefineTrinaryWithVSTScalars(where, OpWhere, Sca) \
231PETE_DefineBinaryWithVSTScalars(Min, FnMin, Sca) \
232PETE_DefineBinaryWithVSTScalars(Max, FnMax, Sca) \
233PETE_DefineBinaryWithVSTScalars(dot, FnDot, Sca) \
234PETE_DefineBinaryWithVSTScalars(dotdot, FnDotDot, Sca) \
235PETE_DefineBinaryWithVSTScalars(outerProduct, FnOuterProduct, Sca) \
236PETE_DefineBinaryWithVSTScalars(cross, FnCross, Sca) \
237PETE_DefineBinaryWithVSTScalars(lt, OpLT, Sca) \
238PETE_DefineBinaryWithVSTScalars(le, OpLE, Sca) \
239PETE_DefineBinaryWithVSTScalars(gt, OpGT, Sca) \
240PETE_DefineBinaryWithVSTScalars(ge, OpGE, Sca) \
241PETE_DefineBinaryWithVSTScalars(eq, OpEQ, Sca) \
242PETE_DefineBinaryWithVSTScalars(ne, OpNE, Sca)
243
247
248#undef PETE_DefineVSTScalar
249
250
251//=========================================================================
252//
253// ASSIGNMENT OPERATIONS
254//
255//=========================================================================
256
257PETE_DefineAssign((a = Min(a,b)),(a = Min(a,b.value)), OpMinAssign)
258PETE_DefineAssign((a = Max(a,b)),(a = Max(a,b.value)), OpMaxAssign)
259PETE_DefineAssign((a = (a&&b)) ,(a = (a&&b.value)) , OpAndAssign)
260PETE_DefineAssign((a = (a||b)) ,(a = (a||b.value)) , OpOrAssign)
261
262
263//=========================================================================
264//
265// MIN and MAX REDUCTIONS
266//
267//=========================================================================
268
269template<class T>
270inline typename T::PETE_Expr_t::PETE_Return_t
271min(const PETE_Expr<T>& expr)
272{
273 typename T::PETE_Expr_t::PETE_Return_t val ;
274 Reduction(val,
276 OpAssign(),
277 OpMinAssign());
278 return val;
279}
280
281template<class T>
282inline typename T::PETE_Expr_t::PETE_Return_t
283max(const PETE_Expr<T>& expr)
284{
285 typename T::PETE_Expr_t::PETE_Return_t val ;
286 Reduction(val,
288 OpAssign(),
289 OpMaxAssign());
290 return val;
291}
292
293//=========================================================================
294//
295// MINMAX REDUCTION
296//
297//=========================================================================
298
299template<class T>
300struct MinMaxHolder {
301 T a;
302 T b;
303 MinMaxHolder() { }
304 MinMaxHolder(const MinMaxHolder<T>& rhs) : a(rhs.a), b(rhs.b) { }
305 const MinMaxHolder<T>& operator=(const MinMaxHolder<T>& rhs) {
306 a = rhs.a;
307 b = rhs.b;
308 return *this;
309 }
310 const MinMaxHolder<T>& operator=(const T& rhs) {
311 T c = rhs;
312 a = c;
313 b = c;
314 return *this;
315 }
316 const MinMaxHolder<T>& operator*=(const MinMaxHolder<T>& rhs) {
317 a = (a < rhs.a ? a : rhs.a);
318 b = (rhs.b < b ? b : rhs.b);
319 return *this;
320 }
321 const MinMaxHolder<T>& operator*=(const T& rhs) {
322 T c = rhs;
323 a = (a < c ? a : c);
324 b = (c < b ? b : c);
325 return *this;
326 }
328 m.put(a);
329 m.put(b);
330 return m;
331 }
333 m.get(a);
334 m.get(b);
335 return m;
336 }
337};
338
339template<class T1, class T2>
340inline void
341minmax(const PETE_Expr<T1>& expr, T2& minval, T2& maxval)
342{
343 typedef typename T1::PETE_Expr_t::PETE_Return_t val_t;
344 MinMaxHolder<val_t> ret;
345 Reduction(ret,
347 OpAssign(),
349 minval = static_cast<T2>(ret.a);
350 maxval = static_cast<T2>(ret.b);
351}
352
353
355//
356// The 'any' function finds if there is any location in the expression
357// where a condition is true.
358//
360
361template<class T, class OP>
363{
364 bool Test;
366 OP Op;
367 AnyHolder() : Test(false), Val(T(0)), Op(OP()) {}
368 AnyHolder(const T& t, OP op) : Test(false), Val(t), Op(op) {}
370 : Test(rhs.Test), Val(rhs.Val), Op(rhs.Op) { }
371 const AnyHolder<T,OP>& operator=(const T& rhs)
372 {
373 if ( PETE_apply(Op,rhs,Val) )
374 Test = true;
375 return *this;
376 }
378 {
379 Test = rhs.Test;
380 Val = rhs.Val;
381 Op = rhs.Op;
382 return *this;
383 }
384 const AnyHolder<T,OP>& operator*=(const T& rhs)
385 {
386 if ( PETE_apply(Op,rhs,Val) )
387 Test = true;
388 return *this;
389 }
391 {
392 Test = (Test || rhs.Test);
393 return *this;
394 }
396 {
397 m.put(Test);
398 return m;
399 }
401 {
402 m.get(Test);
403 return m;
404 }
405};
406
407template<class T1, class T2>
408inline bool
409any(const PETE_Expr<T1>& expr, T2 val)
410{
411 AnyHolder<T2,OpEQ> ret(val,OpEQ());
412 Reduction(ret,
414 OpAssign(),
416 return ret.Test;
417}
418
419template<class T1, class T2, class Op>
420inline bool
421any(const PETE_Expr<T1>& expr, T2 val, Op op)
422{
423 AnyHolder<T2,Op> ret(val,op);
424 Reduction(ret,
426 OpAssign(),
428 return ret.Test;
429}
430
431
432//=========================================================================
433//
434// BOUNDS REDUCTION - find bounding box of Vektor expression
435// for scalars, use minmax
436// for tensors, extend this code to include them as well
437//
438//=========================================================================
439
440template<class T, unsigned int D>
445 BoundsHolder(const BoundsHolder<T,D>& rhs) : a(rhs.a), b(rhs.b) { }
447 a = rhs.a;
448 b = rhs.b;
449 return *this;
450 }
452 Vektor<T,D> c(rhs);
453 a = c;
454 b = c;
455 return *this;
456 }
457 const BoundsHolder<T,D>& operator=(const T& rhs) {
458 Vektor<T,D> c(rhs);
459 a = c;
460 b = c;
461 return *this;
462 }
464 for (unsigned int d=0; d < D; ++d) {
465 a[d] = (a[d] < rhs.a[d] ? a[d] : rhs.a[d]);
466 b[d] = (rhs.b[d] < b[d] ? b[d] : rhs.b[d]);
467 }
468 return *this;
469 }
471 Vektor<T,D> c(rhs);
472 for (unsigned int d=0; d < D; ++d) {
473 a[d] = (a[d] < c[d] ? a[d] : c[d]);
474 b[d] = (c[d] < b[d] ? b[d] : c[d]);
475 }
476 return *this;
477 }
478 const BoundsHolder<T,D>& operator*=(const T& rhs) {
479 Vektor<T,D> c(rhs);
480 for (unsigned int d=0; d < D; ++d) {
481 a[d] = (a[d] < c[d] ? a[d] : c[d]);
482 b[d] = (c[d] < b[d] ? b[d] : c[d]);
483 }
484 return *this;
485 }
487 m.put(a);
488 m.put(b);
489 return m;
490 }
492 m.get(a);
493 m.get(b);
494 return m;
495 }
496};
497
498template<class T1, class T2, unsigned int D>
499inline void
500bounds(const PETE_Expr<T1>& expr, Vektor<T2,D>& minval, Vektor<T2,D>& maxval)
501{
503 Reduction(ret,
505 OpAssign(),
507 minval = ret.a;
508 maxval = ret.b;
509}
510
511
512//=========================================================================
513//
514// OPERATOR()
515//
516//=========================================================================
517
518template<class T, class TP>
521{
522 return a(op.Arg);
523}
524
525
526//=========================================================================
527//
528// MISCELLANEOUS
529//
530//=========================================================================
531
532// When figuring out data dependencies you sometimes need to know
533// if the left hand side of an assignment needs to be read
534// before being it is written.
535// The following trait will set IsAssign to 1 for OpAssign
536// and to zero for all the other assignment functions.
537
538template<class Op> struct OperatorTraits { enum { IsAssign=0 }; };
539template<> struct OperatorTraits<OpAssign> { enum { IsAssign=1 }; };
540
541
543//
544// END OF FILE
545//
547
548#endif // IPPL_EXPRESSIONS_H
549
550/***************************************************************************
551 * $RCSfile: IpplExpressions.h,v $ $Author: adelmann $
552 * $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:28 $
553 * IPPL_VERSION_ID: $Id: IpplExpressions.h,v 1.1.1.1 2003/01/23 07:40:28 adelmann Exp $
554 ***************************************************************************/
FLieGenerator< T, N > imag(const FLieGenerator< std::complex< T >, N > &)
Take imaginary part of a complex generator.
void putMessage(Message &m, const T &t)
Definition: Message.h:549
void getMessage(Message &m, T &t)
Definition: Message.h:572
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
T::PETE_Return_t Reduction(const PETE_Expr< T > &const_expr, CompOp comp_op, AccOp acc_op, NDIndex< D > &loc)
PETE_TBTree< FnDot, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > dot(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
PETE_TBTree< OpGE, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > ge(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
std::complex< double > a
PETEUnaryReturn< T, OpParens< TP > >::type PETE_apply(OpParens< TP > op, const T &a)
FnMax FnDotDot FnCross PETE_DefineBinarySynonym(lt, OpLT) PETE_DefineBinarySynonym(gt
PETE_DefineAssign((a=Min(a, b)),(a=Min(a, b.value)), OpMinAssign) PETE_DefineAssign((a
PETE_TBTree< OpGT, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > gt(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
PETE_TBTree< OpEQ, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > eq(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
short short short short short short int int int int int int long long long long long long float float float float float float double double double double double double OpEQ
#define PETE_DefineScalarMinMax(Sca)
#define PETE_DefineVSTScalar(Sca)
FnArg FnReal FnSign transpose(a))
PETE_TBTree< OpLT, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > lt(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
FnArg real(a))
FnArg FnReal sign(a))
short Max(const short &a, const short &b)
bool any(const PETE_Expr< T1 > &expr, T2 val)
FnMax FnDotDot cross(a, b))
PETE_TBTree< OpNE, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > ne(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
short short short short short short int int int int int int long long long long long long float float float float float float double double double double double double PETE_DefineScalar(std::complex< double >) PETE_DefineBinaryWithScalars(eq
PETE_DefineUnary(Abs,(0< a ? a :-a), FnAbs) inline double PETE_apply(FnAbs
FnMax dotdot(a, b))
short Min(const short &a, const short &b)
#define PETE_DefineIPPLScalar(Sca)
void bounds(const PETE_Expr< T1 > &expr, Vektor< T2, D > &minval, Vektor< T2, D > &maxval)
PETE_DefineBinaryWithScalars(Min, FnMin, short) PETE_DefineBinaryWithScalars(Max
PETE_TBTree< OpLE, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > le(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
void minmax(const PETE_Expr< T1 > &expr, T2 &minval, T2 &maxval)
PETE_DefineBinary(Min,(Min(a, b)), FnMin) PETE_DefineBinary(Max
FnArg FnReal FnSign FnTranspose cofactors(a))
PETE_TBTree< FnOuterProduct, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > outerProduct(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
arg(a))
T det(const AntiSymTenzor< T, 3 > &)
MMatrix< m_complex > complex(MMatrix< double > real)
Definition: MMatrix.cpp:396
m_complex conj(const m_complex &c)
Definition: MVector.h:105
MMatrix< double > & operator*=(MMatrix< double > &m1, MMatrix< double > m2)
Definition: MMatrix.cpp:332
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
BOOST_UBLAS_INLINE V trace(ublas::matrix< V > &e)
Computes the trace of a square matrix.
boost::function< boost::tuple< double, bool >(arguments_t)> type
Definition: function.hpp:21
Definition: Tenzor.h:35
Definition: Vektor.h:32
Message & put(const T &val)
Definition: Message.h:406
Message & get(const T &cval)
Definition: Message.h:476
AnyHolder(const T &t, OP op)
const AnyHolder< T, OP > & operator*=(const AnyHolder< T, OP > &rhs)
const AnyHolder< T, OP > & operator=(const T &rhs)
Message & putMessage(Message &m)
const AnyHolder< T, OP > & operator=(const AnyHolder< T, OP > &rhs)
const AnyHolder< T, OP > & operator*=(const T &rhs)
Message & getMessage(Message &m)
AnyHolder(const AnyHolder< T, OP > &rhs)
const BoundsHolder< T, D > & operator=(const T &rhs)
const BoundsHolder< T, D > & operator*=(const BoundsHolder< T, D > &rhs)
Vektor< T, D > a
Vektor< T, D > b
const BoundsHolder< T, D > & operator=(const BoundsHolder< T, D > &rhs)
Message & getMessage(Message &m)
const BoundsHolder< T, D > & operator*=(const Vektor< T, D > &rhs)
const BoundsHolder< T, D > & operator*=(const T &rhs)
Message & putMessage(Message &m)
BoundsHolder(const BoundsHolder< T, D > &rhs)
const BoundsHolder< T, D > & operator=(const Vektor< T, D > &rhs)
Definition: PETE.h:77
WrappedExpr & PETE_unwrap()
Definition: PETE.h:81