diff options
| -rw-r--r-- | Carpet/CarpetLib/src/vect.cc | 1 | ||||
| -rw-r--r-- | Carpet/CarpetLib/src/vect.hh | 910 | ||||
| -rw-r--r-- | Carpet/CarpetLib/src/vect_helpers.hh | 331 |
3 files changed, 407 insertions, 835 deletions
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc index e42ba04aa..545b0d6d3 100644 --- a/Carpet/CarpetLib/src/vect.cc +++ b/Carpet/CarpetLib/src/vect.cc @@ -18,6 +18,7 @@ void vect<T,D>::input (istream& is) { consume (is, '['); for (int d=0; d<D; ++d) { is >> (*this)[d]; + assert (is.good()); if (d<D-1) { skipws (is); consume (is, ','); diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index 17c75b2b2..29b8a3ba8 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.hh @@ -8,6 +8,8 @@ #include "cctk.h" +#include "vect_helpers.hh" + using namespace std; @@ -109,12 +111,22 @@ public: #if 0 // This creates confusion - /** Constructor from a pointer, i.e.\ a C array. */ + /** Constructor from a pointer, i.e., a C array. */ explicit vect (const T* const x) PURE { for (int d=0; d<D; ++d) elt[d]=x[d]; } #endif +#if 0 + // This leads to an ICE on AIX + template<int E> + operator vect<vect<T,D>,E> () CONST { + vect<vect<T,D>,E> r; + for (int e=0; e<E; ++e) r[e]=*this; + return r; + } +#endif + /** Constructor from a vector with a different type. */ template<typename S> /*explicit*/ vect (const vect<S,D>& a) /*PURE*/ { @@ -225,85 +237,14 @@ public: } // Modifying operators - vect& operator+=(const T& x) { - for (int d=0; d<D; ++d) elt[d]+=x; - return *this; - } - - vect& operator-=(const T& x) { - for (int d=0; d<D; ++d) elt[d]-=x; - return *this; - } - - vect& operator*=(const T& x) { - for (int d=0; d<D; ++d) elt[d]*=x; - return *this; - } - - vect& operator/=(const T& x) { - for (int d=0; d<D; ++d) elt[d]/=x; - return *this; - } - - vect& operator%=(const T& x) { - for (int d=0; d<D; ++d) elt[d]%=x; - return *this; - } - - vect& operator&=(const T& x) { - for (int d=0; d<D; ++d) elt[d]&=x; - return *this; - } - - vect& operator|=(const T& x) { - for (int d=0; d<D; ++d) elt[d]|=x; - return *this; - } - - vect& operator^=(const T& x) { - for (int d=0; d<D; ++d) elt[d]^=x; - return *this; - } - - vect& operator+=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]+=a[d]; - return *this; - } - - vect& operator-=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]-=a[d]; - return *this; - } - - vect& operator*=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]*=a[d]; - return *this; - } - - vect& operator/=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]/=a[d]; - return *this; - } - - vect& operator%=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]%=a[d]; - return *this; - } - - vect& operator&=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]&=a[d]; - return *this; - } - - vect& operator|=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]|=a[d]; - return *this; - } - - vect& operator^=(const vect& a) { - for (int d=0; d<D; ++d) elt[d]^=a[d]; - return *this; - } + DECLARE_MEMBER_OPERATOR_1_REF (operator+=, +=); + DECLARE_MEMBER_OPERATOR_1_REF (operator-=, -=); + DECLARE_MEMBER_OPERATOR_1_REF (operator*=, *=); + DECLARE_MEMBER_OPERATOR_1_REF (operator/=, /=); + DECLARE_MEMBER_OPERATOR_1_REF (operator%=, %=); + DECLARE_MEMBER_OPERATOR_1_REF (operator&=, &=); + DECLARE_MEMBER_OPERATOR_1_REF (operator|=, |=); + DECLARE_MEMBER_OPERATOR_1_REF (operator^=, ^=); // Non-modifying operators @@ -321,189 +262,12 @@ public: return r; } - vect operator+ () const PURE { - vect r; - for (int d=0; d<D; ++d) r[d]=+elt[d]; - return r; - } - - vect operator- () const PURE { - vect r; - for (int d=0; d<D; ++d) r[d]=-elt[d]; - return r; - } - - vect<bool,D> operator! () const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=!elt[d]; - return r; - } - - vect operator~ () const PURE { - vect r; - for (int d=0; d<D; ++d) r[d]=~elt[d]; - return r; - } + DECLARE_MEMBER_OPERATOR_0 (operator+, +) + DECLARE_MEMBER_OPERATOR_0 (operator-, -) + DECLARE_MEMBER_OPERATOR_0 (operator~, ~) + // DECLARE_MEMBER_OPERATOR_0_RET (operator!, !, bool) #if 0 - // These are now defined as template functions - vect operator+ (const T& x) const PURE { - vect r(*this); - r+=x; - return r; - } - - vect operator- (const T& x) const PURE { - vect r(*this); - r-=x; - return r; - } - - vect operator* (const T& x) const PURE { - vect r(*this); - r*=x; - return r; - } - - vect operator/ (const T& x) const PURE { - vect r(*this); - r/=x; - return r; - } - - vect operator% (const T& x) const PURE { - vect r(*this); - r%=x; - return r; - } - - vect operator& (const T& x) const PURE { - vect r(*this); - r&=x; - return r; - } - - vect operator| (const T& x) const PURE { - vect r(*this); - r|=x; - return r; - } - - vect operator^ (const T& x) const PURE { - vect r(*this); - r^=x; - return r; - } - - vect<bool,D> operator&& (const T& x) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]&&x; - return r; - } - - vect<bool,D> operator|| (const T& x) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]||x; - return r; - } - - vect operator+ (const vect& a) const PURE { - vect r(*this); - r+=a; - return r; - } - - vect operator- (const vect& a) const PURE { - vect r(*this); - r-=a; - return r; - } - - vect operator* (const vect& a) const PURE { - vect r(*this); - r*=a; - return r; - } - - vect operator/ (const vect& a) const PURE { - vect r(*this); - r/=a; - return r; - } - - vect operator% (const vect& a) const PURE { - vect r(*this); - r%=a; - return r; - } - - vect operator& (const vect& a) const PURE { - vect r(*this); - r&=a; - return r; - } - - vect operator| (const vect& a) const PURE { - vect r(*this); - r|=a; - return r; - } - - vect operator^ (const vect& a) const PURE { - vect r(*this); - r^=a; - return r; - } - - vect<bool,D> operator&& (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]&&a[d]; - return r; - } - - vect<bool,D> operator|| (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]||a[d]; - return r; - } - - vect<bool,D> operator== (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]==a[d]; - return r; - } - - vect<bool,D> operator!= (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]!=a[d]; - return r; - } - - vect<bool,D> operator< (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]<a[d]; - return r; - } - - vect<bool,D> operator<= (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]<=a[d]; - return r; - } - - vect<bool,D> operator> (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]>a[d]; - return r; - } - - vect<bool,D> operator>= (const vect& a) const PURE { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=elt[d]>=a[d]; - return r; - } -#endif - /** This corresponds to the ?: operator. Return a vector with the elements set to either a[i] or b[i], depending on whether (*this)[i] is true or not. */ @@ -513,6 +277,7 @@ public: for (int d=0; d<D; ++d) r[d]=elt[d]?a[d]:b[d]; return r; } +#endif // Iterators #if 0 @@ -572,66 +337,6 @@ inline vect<vect<T,D>,DD> xpose (vect<vect<T,DD>,D> const & a) { return r; } -/** Return the element-wise absolute value. */ -template<typename T,int D> -inline vect<T,D> abs (const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> abs (const vect<T,D>& a) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=abs(a[d]); - return r; -} - -/** Return the element-wise ceiling. */ -template<typename T,int D> -inline vect<T,D> ceil (const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> ceil (const vect<T,D>& a) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=ceil(a[d]); - return r; -} - -/** Return the element-wise floor. */ -template<typename T,int D> -inline vect<T,D> floor (const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> floor (const vect<T,D>& a) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=floor(a[d]); - return r; -} - -/** Return the element-wise maximum of two vectors. */ -template<typename T,int D> -inline vect<T,D> max (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> max (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=max(a[d],b[d]); - return r; -} - -/** Return the element-wise minimum of two vectors. */ -template<typename T,int D> -inline vect<T,D> min (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> min (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=min(a[d],b[d]); - return r; -} - -/** Return the element-wise power of two vectors. */ -template<typename T,typename U,int D> -inline vect<T,D> pow (const vect<T,D>& a, const vect<U,D>& b) PURE; -template<typename T,typename U,int D> -inline vect<T,D> pow (const vect<T,D>& a, const vect<U,D>& b) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=pow(a[d],b[d]); - return r; -} - /** Return the element-wise integer power of two vectors. */ template<typename T,int D> inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b) PURE; @@ -642,445 +347,51 @@ inline vect<T,D> ipow (const vect<T,D>& a, const vect<int,D>& b) { return r; } -template<typename T,int D> -inline vect<T,D> operator+ (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator+ (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r+=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator+ (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator+ (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r+=a; - return r; -} +DECLARE_FUNCTION_1 (abs) +DECLARE_FUNCTION_1 (ceil) +DECLARE_FUNCTION_1 (floor) -template<typename T,int D> -inline vect<T,D> operator- (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator- (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r-=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator- (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator- (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r-=a; - return r; -} +DECLARE_OPERATOR_1_RET (operator!, !, bool) -template<typename T,int D> -inline vect<T,D> operator* (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator* (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r*=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator* (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator* (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r*=a; - return r; -} +DECLARE_FUNCTION_2 (max) +DECLARE_FUNCTION_2 (min) +DECLARE_FUNCTION_2 (pow) -template<typename T,int D> -inline vect<T,D> operator/ (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator/ (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r/=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator/ (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator/ (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r/=a; - return r; -} +DECLARE_OPERATOR_2 (operator+, +) +DECLARE_OPERATOR_2 (operator-, -) +DECLARE_OPERATOR_2 (operator*, *) +DECLARE_OPERATOR_2 (operator/, /) +DECLARE_OPERATOR_2 (operator%, %) +DECLARE_OPERATOR_2 (operator&, &) +DECLARE_OPERATOR_2 (operator|, |) +DECLARE_OPERATOR_2 (operator^, ^) -template<typename T,int D> -inline vect<T,D> operator% (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator% (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r%=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator% (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator% (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r%=a; - return r; -} +DECLARE_OPERATOR_2_RET (operator&&, &&, bool) +DECLARE_OPERATOR_2_RET (operator||, ||, bool) +DECLARE_OPERATOR_2_RET (operator==, ==, bool) +DECLARE_OPERATOR_2_RET (operator!=, !=, bool) +DECLARE_OPERATOR_2_RET (operator< , < , bool) +DECLARE_OPERATOR_2_RET (operator<=, <=, bool) +DECLARE_OPERATOR_2_RET (operator> , > , bool) +DECLARE_OPERATOR_2_RET (operator>=, >=, bool) -template<typename T,int D> -inline vect<T,D> operator| (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator| (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r|=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator| (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator| (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r|=a; - return r; -} -template<typename T,int D> -inline vect<T,D> operator& (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator& (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r&=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator& (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator& (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r&=a; - return r; -} -template<typename T,int D> -inline vect<T,D> operator^ (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator^ (const vect<T,D>& a, const T& x) { - vect<T,D> r(a); - r^=x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator^ (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator^ (const T& x, const vect<T,D>& a) { - vect<T,D> r(x); - r^=a; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator&& (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator&& (const vect<T,D>& a, const T& x) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]&&x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator&& (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator&& (const T& x, const vect<T,D>& a) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=x&&a[d]; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator|| (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<T,D> operator|| (const vect<T,D>& a, const T& x) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]||x; - return r; -} -template<typename T,int D> -inline vect<T,D> operator|| (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<T,D> operator|| (const T& x, const vect<T,D>& a) { - vect<T,D> r; - for (int d=0; d<D; ++d) r[d]=x||a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator== (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator== (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]==x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator== (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator== (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x==a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator!= (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator!= (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]!=x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator!= (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator!= (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x!=a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator< (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator< (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]<x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator< (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator< (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x<a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator<= (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator<= (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]<=x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator<= (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator<= (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x<=a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator> (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator> (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]>x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator> (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator> (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x>a[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator>= (const vect<T,D>& a, const T& x) PURE; -template<typename T,int D> -inline vect<bool,D> operator>= (const vect<T,D>& a, const T& x) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]>=x; - return r; -} -template<typename T,int D> -inline vect<bool,D> operator>= (const T& x, const vect<T,D>& a) PURE; -template<typename T,int D> -inline vect<bool,D> operator>= (const T& x, const vect<T,D>& a) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=x>=a[d]; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator+ (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator+ (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r+=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator- (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator- (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r-=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator* (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator* (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r*=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator/ (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator/ (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r/=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator% (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator% (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r%=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator& (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator& (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r&=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator| (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator| (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r|=b; - return r; -} - -template<typename T,int D> -inline vect<T,D> operator^ (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<T,D> operator^ (const vect<T,D>& a, const vect<T,D>& b) { - vect<T,D> r(a); - r^=b; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator&& (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator&& (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]&&b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator|| (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator|| (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]||b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator== (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator== (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]==b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator!= (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator!= (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]!=b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator< (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator< (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]<b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator<= (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator<= (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]<=b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator> (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator> (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]>b[d]; - return r; -} - -template<typename T,int D> -inline vect<bool,D> operator>= (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline vect<bool,D> operator>= (const vect<T,D>& a, const vect<T,D>& b) { - vect<bool,D> r; - for (int d=0; d<D; ++d) r[d]=a[d]>=b[d]; - return r; -} +// Reduction operators +// Identity +#define id(x) (x) +DECLARE_REDUCTION_OPERATOR_1_T_RET (all,true ,&=,id,bool,bool) +DECLARE_REDUCTION_OPERATOR_1_T_RET (any,false,|=,id,bool,bool) -// Reduction operators +DECLARE_REDUCTION_FUNCTION_1 (maxval,a[0],max,id) +DECLARE_REDUCTION_FUNCTION_1 (minval,a[0],min,id) +DECLARE_REDUCTION_OPERATOR_1 (prod,1,*=,id) +DECLARE_REDUCTION_OPERATOR_1 (sum,0,+=,id) -/** Return true iff any of the elements are true (boolean sum). */ -template<int D> -inline bool any (const vect<bool,D>& a) PURE; -template<int D> -inline bool any (const vect<bool,D>& a) { - bool r(false); - for (int d=0; d<D; ++d) r|=a[d]; - return r; -} - -/** Return true iff all of the elements are true (boolean product). */ -template<int D> -inline bool all (const vect<bool,D>& a) PURE; -template<int D> -inline bool all (const vect<bool,D>& a) { - bool r(true); - for (int d=0; d<D; ++d) r&=a[d]; - return r; -} +DECLARE_REDUCTION_OPERATOR_2 (dot ,0,+=,*,id ) +DECLARE_REDUCTION_OPERATOR_2 (hypot,0,+=,*,sqrt) /** Count the number of elements in the vector. */ template<typename T,int D> @@ -1090,44 +401,12 @@ inline int count (const vect<T,D>& a) { return D; } -/** Return the dot product of two vectors. */ -template<typename T,int D> -inline T dot (const vect<T,D>& a, const vect<T,D>& b) PURE; -template<typename T,int D> -inline T dot (const vect<T,D>& a, const vect<T,D>& b) { - T r(0); - for (int d=0; d<D; ++d) r+=a[d]*b[d]; - return r; -} - -/** Return the Euklidean length. */ -template<typename T,int D> -inline T hypot (const vect<T,D>& a) PURE; -template<typename T,int D> -inline T hypot (const vect<T,D>& a) { - return sqrt(dot(a,a)); -} - -/** Return the maximum element. */ -template<typename T,int D> -inline T maxval (const vect<T,D>& a) PURE; -template<typename T,int D> -inline T maxval (const vect<T,D>& a) { - assert (D>0); - T r(a[0]); - for (int d=1; d<D; ++d) r=max(r,a[d]); - return r; -} - -/** Return the minimum element. */ +/** Return the size (number of elements) of the vector. */ template<typename T,int D> -inline T minval (const vect<T,D>& a) PURE; +inline int size (const vect<T,D>& a) CONST; template<typename T,int D> -inline T minval (const vect<T,D>& a) { - assert (D>0); - T r(a[0]); - for (int d=1; d<D; ++d) r=min(r,a[d]); - return r; +inline int size (const vect<T,D>& a) { + return D; } /** Return the index of the first maximum element. */ @@ -1152,33 +431,7 @@ inline int minloc (const vect<T,D>& a) { return r; } -/** Return the product of the elements. */ -template<typename T,int D> -inline T prod (const vect<T,D>& a) PURE; -template<typename T,int D> -inline T prod (const vect<T,D>& a) { - T r(1); - for (int d=0; d<D; ++d) r*=a[d]; - return r; -} -/** Return the size (number of elements) of the vector. */ -template<typename T,int D> -inline int size (const vect<T,D>& a) CONST; -template<typename T,int D> -inline int size (const vect<T,D>& a) { - return D; -} - -/** Return the sum of the elements. */ -template<typename T,int D> -inline T sum (const vect<T,D>& a) PURE; -template<typename T,int D> -inline T sum (const vect<T,D>& a) { - T r(0); - for (int d=0; d<D; ++d) r+=a[d]; - return r; -} // Higher order functions @@ -1308,30 +561,6 @@ vect (const T& x, const T& y, const T& z, const T& t) { //////////////////////////////////////////////////////////////////////////////// -// Special versions for vect<vect> - -template<typename T,int D,int E> -inline vect<vect<bool,D>,E> operator== (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) PURE; -template<typename T,int D,int E> -inline vect<vect<bool,D>,E> operator== (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) { - vect<vect<bool,D>,E> r; - for (int e=0; e<E; ++e) r[e]=a[e]==b[e]; - return r; -} - -/** Return true iff all of the elements are true (boolean product). */ -template<int D,int E> -inline vect<bool,E> all (const vect<vect<bool,D>,E>& a) PURE; -template<int D,int E> -inline vect<bool,E> all (const vect<vect<bool,D>,E>& a) { - vect<bool,E> r; - for (int e=0; e<E; ++e) r[e]=all(a[e]); - return r; -} - - -//////////////////////////////////////////////////////////////////////////////// - // Specialise some constructors for lower dimensions @@ -1365,6 +594,17 @@ inline vect<CCTK_REAL,3>& vect<CCTK_REAL,3>::operator%=(const vect<CCTK_REAL,3>& return *this; } +template<> +inline vect<CCTK_REAL,3> operator%(const vect<CCTK_REAL,3>& a, const vect<CCTK_REAL,3>& b) { + vect<CCTK_REAL,3> r; + for (int d=0; d<3; ++d) { + r[d]=fmod(a[d],b[d]); + if (r[d]>b[d]*(CCTK_REAL)(1.0-1.0e-10)) r[d]=(CCTK_REAL)0; + if (r[d]<b[d]*(CCTK_REAL)( 1.0e-10)) r[d]=(CCTK_REAL)0; + } + return r; +} + #endif // VECT_HH diff --git a/Carpet/CarpetLib/src/vect_helpers.hh b/Carpet/CarpetLib/src/vect_helpers.hh new file mode 100644 index 000000000..afcce6e7e --- /dev/null +++ b/Carpet/CarpetLib/src/vect_helpers.hh @@ -0,0 +1,331 @@ +#ifndef VECT_HELPERS_HH +#define VECT_HELPERS_HH + + + +// Declare a member operator which takes 0 arguments + +#define DECLARE_MEMBER_OPERATOR_0(fn,op) \ + \ + vect fn () const PURE { \ + vect r; \ + for (int d=0; d<D; ++d) r[d]=op elt[d]; \ + return r; \ + } + + + +// Declare a member operator which takes 0 arguments and returns type +// R + +#define DECLARE_MEMBER_OPERATOR_0_RET(fn,op,R) \ + \ + vect<R,D> fn () const PURE { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=op elt[d]; \ + return r; \ + } + + + +// Declare a member operator which takes 1 argument and returns a +// reference + +#define DECLARE_MEMBER_OPERATOR_1_REF(fn,op) \ + \ + vect& fn (const T& x) { \ + for (int d=0; d<D; ++d) elt[d] op x; \ + return *this; \ + } \ + \ + vect& fn (const vect& a) { \ + for (int d=0; d<D; ++d) elt[d] op a[d]; \ + return *this; \ + } + + + +// Declare a function which takes 1 argument and returns type R + +#define DECLARE_FUNCTION_1_RET(fn,R) \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=fn(a[d]); \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<R,D> fn (const vect<vect<T,D>,E>& a) PURE; \ + template<typename T,int D,int E> \ + inline vect<R,D> fn (const vect<vect<T,D>,E>& a) \ + { \ + vect<R,D> r; \ + for (int e=0; e<E; ++e) r[e]=fn(a[e]); \ + return r; \ + } + + + +// Declare a function which takes 1 argument + +#define DECLARE_FUNCTION_1(fn) DECLARE_FUNCTION_1_RET(fn,T) + + + +// Declare a function which takes 2 arguments and returns type R + +#define DECLARE_FUNCTION_2_RET(fn,R) \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const vect<T,D>& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const vect<T,D>& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=fn(a[d],b[d]); \ + return r; \ + } \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a, const vect<T,D>& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a, const vect<T,D>& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=fn(a,b[d]); \ + return r; \ + } \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const T& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const T& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=fn(a[d],b); \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=fn(a[e],b[e]); \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a, const vect<vect<T,D>,E>& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a, const vect<vect<T,D>,E>& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=fn(a,b[e]); \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const T& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const T& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=fn(a[e],b); \ + return r; \ + } + + + +// Declare a function which takes 2 arguments + +#define DECLARE_FUNCTION_2(fn) DECLARE_FUNCTION_2_RET(fn,T) + + + +// Declare an operator which takes 1 argument and returns type R + +#define DECLARE_OPERATOR_1_RET(fn,op,R) \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=op a[d]; \ + return r; \ + } \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=op a; \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=op a[e]; \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=op a; \ + return r; \ + } + + + +// Declare an operator which takes 2 arguments and returns type R + +#define DECLARE_OPERATOR_2_RET(fn,op,R) \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const vect<T,D>& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const vect<T,D>& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=a[d] op b[d]; \ + return r; \ + } \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a, const vect<T,D>& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const T& a, const vect<T,D>& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=a op b[d]; \ + return r; \ + } \ + \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const T& b) PURE; \ + template<typename T,int D> \ + inline vect<R,D> fn (const vect<T,D>& a, const T& b) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=a[d] op b; \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const vect<vect<T,D>,E>& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=a[e] op b[e]; \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a, const vect<vect<T,D>,E>& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const T& a, const vect<vect<T,D>,E>& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=a op b[e]; \ + return r; \ + } \ + \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const T& b) PURE; \ + template<typename T,int D,int E> \ + inline vect<vect<R,D>,E> fn (const vect<vect<T,D>,E>& a, const T& b) \ + { \ + vect<vect<R,D>,E> r; \ + for (int e=0; e<E; ++e) r[e]=a[e] op b; \ + return r; \ + } + + + +// Declare an operator which takes 2 arguments + +#define DECLARE_OPERATOR_2(fn,op) DECLARE_OPERATOR_2_RET(fn,op,T) + + + +// Declare a reduction function which takes 1 argument of type T and +// returns type R + +#define DECLARE_REDUCTION_OPERATOR_1_T_RET(fn,init,op,final,T,R) \ + \ + template<typename U,int D> \ + inline vect<R,D> fn (const vect<U,D>& a) PURE; \ + template<typename U,int D> \ + inline vect<R,D> fn (const vect<U,D>& a) { \ + vect<R,D> r; \ + for (int d=0; d<D; ++d) r[d]=fn(a[d]); \ + return r; \ + } \ + \ + template<int D> \ + inline R fn (const vect<T,D>& a) PURE; \ + template<int D> \ + inline R fn (const vect<T,D>& a) { \ + R r(init); \ + for (int d=0; d<D; ++d) r op a[d]; \ + return final(r); \ + } + + + +// Declare a reduction function which takes 1 argument + +#define DECLARE_REDUCTION_OPERATOR_1(fn,init,op,final) \ + \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a) PURE; \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a) { \ + T r(init); \ + for (int d=0; d<D; ++d) r op a[d]; \ + return final(r); \ + } + + + +// Declare a reduction function which takes 1 argument + +#define DECLARE_REDUCTION_FUNCTION_1(fn,init,op,final) \ + \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a) PURE; \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a) { \ + T r(init); \ + for (int d=0; d<D; ++d) op(r,a[d]); \ + return final(r); \ + } + + + +// Declare a reduction function which takes 2 arguments + +#define DECLARE_REDUCTION_OPERATOR_2(fn,init,op,op2,final) \ + \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a, const vect<T,D>& b) PURE; \ + template<typename T,int D> \ + inline T fn (const vect<T,D>& a, const vect<T,D>& b) { \ + T r(init); \ + for (int d=0; d<D; ++d) r op (a[d] op2 b[d]); \ + return final(r); \ + } + + + +#endif // #ifndef VECT_HELPERS_HH |