Vektor.h

Go to the documentation of this file.

// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
 *
 *
 * Visit http://people.web.psi.ch/adelmann/ for more details
 *
 ***************************************************************************/
 
#ifndef VEKTOR_H
#define VEKTOR_H
 
// include files
#include "Utility/PAssert.h"
#include "Message/Message.h"
#include "PETE/IpplExpressions.h"
#include "AppTypes/TSVMeta.h"
 
#include <cmath>
#include <iostream>
#include <iomanip>
 
//
// Definition of class Vektor.
//
 
template<class T, unsigned D>
class Vektor
{
public:
 
  typedef T Element_t;
  enum { ElemDim = 1 };
  enum { Size = D };
 
  // Default Constructor initializes to zero.
  Vektor() {
    TSV_MetaAssignScalar<Vektor<T,D>,T,OpAssign>::apply(*this,T(0));
  }
 
  // Copy Constructor
  Vektor(const Vektor<T,D> &rhs) {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T,D> ,OpAssign >::apply(*this,rhs);
  }
 
  // Templated Vektor constructor.
  template<class T1, unsigned D1>
  Vektor(const Vektor<T1,D1> &rhs) {
    for (unsigned d=0; d<D; ++d)
      X[d] = (d < D1) ? rhs[d] : T1(0);
  }
 
  // Constructor from a single T
  Vektor(const T& x00) {
    TSV_MetaAssignScalar<Vektor<T,D>,T,OpAssign>::apply(*this,x00);
  }
 
  Vektor(std::initializer_list<T> init) {
        unsigned i = 0;
        for(auto it = init.begin(); it != init.end() && i < D; ++it, ++i) {
            X[i] = *it; // copy at most D elements to avoid out-of-bounds access
        }
        for(; i < D; ++i) {
            X[i] = T(0); // zero-fill remaining components if initializer list is shorter
        }
    }
 
  // Destructor
  ~Vektor() { }
 
  // Assignment Operators
  const Vektor<T,D>& operator=(const Vektor<T,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T,D> ,OpAssign> :: apply(*this,rhs);
    return *this;
  }
  template<class T1>
  const Vektor<T,D>& operator=(const Vektor<T1,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T1,D> ,OpAssign> :: apply(*this,rhs);
    return *this;
  }
  const Vektor<T,D>& operator=(const T& rhs)
  {
    TSV_MetaAssignScalar< Vektor<T,D> , T ,OpAssign > :: apply(*this,rhs);
    return *this;
  }
 
  // Accumulation Operators
  template<class T1>
  Vektor<T,D>& operator+=(const Vektor<T1,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T1,D> , OpAddAssign > :: apply(*this,rhs);
    return *this;
  }
  Vektor<T,D>& operator+=(const T& rhs)
  {
    TSV_MetaAssignScalar< Vektor<T,D> , T , OpAddAssign > :: apply(*this,rhs);
    return *this;
  }
 
  template<class T1>
  Vektor<T,D>& operator-=(const Vektor<T1,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T1,D> , OpSubtractAssign > :: apply(*this,rhs);
    return *this;
  }
  Vektor<T,D>& operator-=(const T& rhs)
  {
    TSV_MetaAssignScalar< Vektor<T,D> , T , OpSubtractAssign > :: apply(*this,rhs);
    return *this;
  }
 
  template<class T1>
  Vektor<T,D>& operator*=(const Vektor<T1,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T1,D> , OpMultipplyAssign > :: apply(*this,rhs);
    return *this;
  }
  Vektor<T,D>& operator*=(const T& rhs)
  {
    TSV_MetaAssignScalar< Vektor<T,D> , T , OpMultipplyAssign > :: apply(*this,rhs);
    return *this;
  }
 
  template<class T1>
  Vektor<T,D>& operator/=(const Vektor<T1,D> &rhs)
  {
    TSV_MetaAssign< Vektor<T,D> , Vektor<T1,D> , OpDivideAssign > ::
      apply(*this,rhs);
    return *this;
  }
  Vektor<T,D>& operator/=(const T& rhs)
  {
    TSV_MetaAssignScalar< Vektor<T,D> , T , OpDivideAssign > ::
      apply(*this,rhs);
    return *this;
  }
 
  // Get and Set Operations
  Element_t& operator[](unsigned int i);
 
  Element_t operator[](unsigned int i) const;
 
  Element_t& operator()(unsigned int i);
 
  Element_t operator()( unsigned int i) const;
 
  // Comparison operators.
  bool operator==(const Vektor<T,D>& that) const {
    return TSV_MetaCompareArrays<T,T,D>::apply(X,that.X);
  }
  bool operator!=(const Vektor<T,D>& that) const {
    return !(*this == that);
  }
 
  //----------------------------------------------------------------------
  // parallel communication
  Message& putMessage(Message& m) const {
    m.setCopy(true);
    ::putMessage(m, X, X + D);
    return m;
  }
 
  Message& getMessage(Message& m) {
    ::getMessage(m, X, X + D);
    return m;
  }
 
private:
 
  // Just store D elements of type T.
    T X[D];
 
};
 
template<class T, unsigned D>
typename Vektor<T,D>::Element_t& Vektor<T,D>::operator[](unsigned int i)
{
  PAssert (i<D);
  if (i >= D) i = D-1; // clamp index to avoid undefined behavior when assertions are disabled
  return X[i];
}
 
template<class T, unsigned D>
typename Vektor<T,D>::Element_t Vektor<T,D>::operator[](unsigned int i) const
{
  PAssert (i<D);
  if (i >= D) i = D-1; // defensive fallback for release builds without runtime bounds checking
  return X[i];
}
 
template<class T, unsigned D>
typename Vektor<T,D>::Element_t& Vektor<T,D>::operator()(unsigned int i)
{
  PAssert (i<D);
  return X[i];
}
 
template<class T, unsigned D>
typename Vektor<T,D>::Element_t Vektor<T,D>::operator()( unsigned int i) const
{
  PAssert (i<D);
  return X[i];
}
 
//
// Unary Operators
//
 
//----------------------------------------------------------------------
// unary operator-
template<class T, unsigned D>
inline Vektor<T,D> operator-(const Vektor<T,D> &op)
{
  return TSV_MetaUnary< Vektor<T,D> , OpUnaryMinus > :: apply(op);
}
 
//----------------------------------------------------------------------
// unary operator+
template<class T, unsigned D>
inline const Vektor<T,D> &operator+(const Vektor<T,D> &op)
{
  return op;
}
 
//
// Binary Operators
//
 
//
// Elementwise operators.
//
 
TSV_ELEMENTWISE_OPERATOR(Vektor,operator+,OpAdd)
TSV_ELEMENTWISE_OPERATOR(Vektor,operator-,OpSubtract)
TSV_ELEMENTWISE_OPERATOR(Vektor,operator*,OpMultipply)
TSV_ELEMENTWISE_OPERATOR(Vektor,operator/,OpDivide)
TSV_ELEMENTWISE_OPERATOR(Vektor,Min,FnMin)
TSV_ELEMENTWISE_OPERATOR(Vektor,Max,FnMax)
 
//----------------------------------------------------------------------
// dot product
//----------------------------------------------------------------------
 
template < class T1, class T2, unsigned D >
inline typename PETEBinaryReturn<T1,T2,OpMultipply>::type
dot(const Vektor<T1,D> &lhs, const Vektor<T2,D> &rhs)
{
  return TSV_MetaDot< Vektor<T1,D> , Vektor<T2,D> > :: apply(lhs,rhs);
}
 
//----------------------------------------------------------------------
// cross product
//----------------------------------------------------------------------
 
template < class T1, class T2, unsigned D >
inline Vektor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
cross(const Vektor<T1,D> &lhs, const Vektor<T2,D> &rhs)
{
  return TSV_MetaCross< Vektor<T1,D> , Vektor<T2,D> > :: apply(lhs,rhs);
}
 
//----------------------------------------------------------------------
// euclidean norm
//----------------------------------------------------------------------
template < class T, unsigned D >
inline double
euclidean_norm(const Vektor<T, D> &a)
{
  return std::sqrt(dot(a, a));
}
 
 
//----------------------------------------------------------------------
// I/O
template<class T, unsigned D>
inline std::ostream& operator<<(std::ostream& out, const Vektor<T,D>& rhs)
{
  std::streamsize sw = out.width();
  out << std::setw(1);
  if (D >= 1) {
    out << "( ";
    for (unsigned int i=0; i<D - 1; i++)
      out << std::setw(sw) << rhs[i] << " , ";
    out << std::setw(sw) << rhs[D - 1] << " )";
  } else {
    out << "( " << std::setw(sw) << rhs[0] << " )";
  }
 
  return out;
}
 
#endif // VEKTOR_H