Vec.h

#ifndef PSCF_VEC_H
#define PSCF_VEC_H
 
/*
* PSCF - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <util/global.h>
#include <util/accumulators/setToZero.h>
#include <iostream>
 
using namespace Util;
 
namespace Pscf
{
 
   template <int D, typename T>
   class Vec
   {
 
   public:
 
 
 
      Vec<D, T>();
 
      Vec<D, T>(const Vec<D, T>& v);
 
      explicit Vec<D, T>(T const *  v);
 
      explicit Vec<D, T>(T s);
 
 
 
 
      Vec<D, T>& operator = (const Vec<D, T>& v);
 
      Vec<D, T>& operator = (T s);
 
      Vec<D, T>& setToZero();
 
 
 
      void operator += (const Vec<D, T>& dv);
 
      void operator -= (const Vec<D, T>& dv);
 
      void operator += (T s);
 
      void operator -= (T s);
 
      void operator *= (T s);
 
 
 
      const T& operator [] (int i) const;
 
      T& operator [] (int i);
 
 
 
      Vec<D, T>& add(const Vec<D, T>& v1, const Vec<D, T>& v2);
 
      Vec<D, T>& subtract(const Vec<D, T>& v1, const Vec<D, T>& v2);
 
      Vec<D, T>& multiply(const Vec<D, T>& v, T s);
 
      Vec<D, T>& negate(const Vec<D, T>& v);
 
      Vec<D, T>& negate();
 
 
      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);
 
   private:
 
      static const int Width = 25;
 
      static const int Precision = 17;
 
      T elem_[D];
 
   };
 
   // Associated functions
 
   template <int D, typename T>
   inline 
   T dot(Vec<D, T> const & v1, Vec<D, T> const & v2)
   {
      T value;
      setToZero(value);
      for (int i = 0; i < D; ++i) {
         value += v1[i]*v2[i];
      }
      return value;
   }
 
   template <int D, typename T>
   inline 
   Vec<D, T> operator + (Vec<D, T> const & v1, Vec<D, T> const & v2)
   {
      Vec<D, T> value;
      value.add(v1, v2);
      return value;
   }
 
   // Inline method definitions
 
   /*
   * Default constructor
   */
   template <int D, typename T> 
   inline Vec<D, T>::Vec()
   {}
 
   /*
   * Copy constructor
   */
   template <int D, typename T> 
   inline Vec<D, T>::Vec(const Vec<D, T>& v)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v.elem_[i];
      }
   }
 
   /*
   * Constructor, from C-Array.
   */
   template <int D, typename T> 
   inline Vec<D, T>::Vec(T const * v)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v[i];
      }
   }
 
   /*
   * Constructor, initialize all elements to a scalar value s.
   */
   template <int D, typename T> 
   inline Vec<D, T>::Vec(T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = s;
      }
   }
 
   /*
   * Assignment.
   */
   template <int D, typename T> 
   inline Vec<D, T>& Vec<D, T>::operator = (const Vec<D, T>& v)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v.elem_[i];
      }
      return *this;
   }
 
   /*
   * Assign all elements to a common scalar value.
   */
   template <int D, typename T> 
   inline Vec<D, T>& Vec<D, T>::operator = (T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = s;
      }
      return *this;
   }
 
   /*
   * Set all elements of a 3D vector to zero.
   */
   template <int D, typename T> 
   inline Vec<D, T>& Vec<D, T>::setToZero()
   {
      for (int i = 0; i < D; ++i) {
         setToZero(elem_[i]); 
      }
      return *this;
   }
 
   /*
   * Add vector dv to this vector.
   */
   template <int D, typename T> 
   inline void Vec<D, T>::operator += (const Vec<D, T>& dv)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] += dv.elem_[i];
      }
   }
 
   /*
   * Subtract vector dv from this vector.
   */
   template <int D, typename T> 
   inline void Vec<D, T>::operator -= (const Vec<D, T>& dv)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] -= dv.elem_[i];
      }
   }
 
   /*
   * Add a common scalar to all components.
   */
   template <int D, typename T> 
   inline void Vec<D, T>::operator += (T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] += s;
      }
   }
 
   /*
   * Subtract a common scalar from all components.
   */
   template <int D, typename T> 
   inline void Vec<D, T>::operator -= (T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] -= s;
      }
   }
 
   /*
   * Multiply this vector by scalar s.
   */
   template <int D, typename T> 
   inline void Vec<D, T>::operator *= (T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] *= s;
      }
   }
 
   /*
   * Return one Cartesian element by value.
   */
   template <int D, typename T>
   inline const T& Vec<D, T>::operator [] (int i) const
   { 
      assert(i >=0); 
      assert(i < D); 
      return elem_[i]; 
   }
 
   /*
   * Return a reference to one element of the vector.
   */
   template <int D, typename T>
   inline T& Vec<D, T>::operator [] (int i)
   {
      assert(i >=0); 
      assert(i < D); 
      return elem_[i]; 
   }
 
   /*
   * Add vectors v1 and v2.
   *
   * Upon return, *this = v1 + v2.
   */
   template <int D, typename T>
   inline
   Vec<D, T>& Vec<D, T>::add(Vec<D, T> const & v1, 
                             Vec<D, T> const & v2)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v1.elem_[i] + v2.elem_[i];
      }
      return *this;
   }
 
   /*
   * Subtract vector v2 from v1.
   *
   * Upon return, *this = v1 - v2.
   * \return modified invoking vector
   */
   template <int D, typename T>
   inline
   Vec<D, T>& Vec<D, T>::subtract(Vec<D, T> const & v1, Vec<D, T> const & v2)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v1.elem_[i] - v2.elem_[i];
      }
      return *this;
   }
 
   /*
   * Multiply a vector v by a scalar s.
   *
   * Upon return, *this = v*s.
   */
   template <int D, typename T>
   inline 
   Vec<D, T>& Vec<D, T>::multiply(Vec<D, T> const & v, T s)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = v.elem_[i]*s;
      }
      return *this;
   }
 
   /*
   * Compute and return negation of a vector.
   *
   * Upon return, *this = -v.
   */
   template <int D, typename T>
   inline 
   Vec<D, T>& Vec<D, T>::negate(Vec<D, T> const & v)
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = -v.elem_[i];
      }
      return *this;
   }
 
   /*
   * Negate (reverse sign) of this vector.
   *
   * Upon return, *this = -v.
   */
   template <int D, typename T>
   inline 
   Vec<D, T>& Vec<D, T>::negate()
   {
      for (int i = 0; i < D; ++i) {
         elem_[i] = -elem_[i];
      }
      return *this;
   }
 
   /*
   * Serialize to/from an archive.
   */
   template <int D, typename T>
   template <class Archive>
   inline 
   void Vec<D, T>::serialize(Archive& ar, const unsigned int version)
   { 
      for (int i = 0; i < D; ++i) {
         ar & elem_[i];
      }
   }
 
}
#endif