Typed vector and tensor

regtest/basic/rt-make-0/config

@@ -1 +1,2 @@
 type=make
+#should we rename this to something meaningful? like "rt-VectorTensor"?

regtest/basic/rt-make-4/main.cpp

@@ -50,7 +50,9 @@ int main () {
     PLMD::TensorGeneric<4,4> mat;
     PLMD::TensorGeneric<1,4> evec;
     PLMD::VectorGeneric<4> eval_underlying;
-
+
+    //The both of the two following lines are scary, but on my machine are equivalent
+    //PLMD::Vector1d* eval=reinterpret_cast<PLMD::Vector1d*>(eval_underlying.data());
     auto eval = new(&eval_underlying[0]) PLMD::VectorGeneric<1>;
 
     mat[1][0]=mat[0][1]=3.0;

src/tools/Communicator.h

@@ -73,13 +73,13 @@ class Communicator {
 /// Init from reference
     template <typename T> explicit Data(T&p): pointer(&p), size(1), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from pointer to VectorGeneric
-    template <unsigned n> explicit Data(VectorGeneric<n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T, unsigned n> explicit Data(VectorTyped<T,n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from reference to VectorGeneric
-    template <unsigned n> explicit Data(VectorGeneric<n> &p): pointer(&p), size(n), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-/// Init from pointer to TensorGeneric
-    template <unsigned n,unsigned m> explicit Data(TensorGeneric<n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-/// Init from reference to TensorGeneric
-    template <unsigned n,unsigned m> explicit Data(TensorGeneric<n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T, unsigned n> explicit Data(VectorTyped<T,n> &p): pointer(&p), size(n), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+/// Init from pointer to TensorTyped
+    template <typename T, unsigned n,unsigned m> explicit Data(TensorTyped<T,n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+/// Init from reference to TensorTyped
+    template <typename T, unsigned n,unsigned m> explicit Data(TensorTyped<T,n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from reference to std::vector
     template <typename T> explicit Data(std::vector<T>&v) {
       Data d(v.data(),v.size());
@@ -121,10 +121,10 @@ class Communicator {
     MPI_Datatype type;
     template <typename T> explicit ConstData(const T*p,int s): pointer(p), size(s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
     template <typename T> explicit ConstData(const T&p): pointer(&p), size(1), nbytes(sizeof(T)), type(getMPIType<T>()) {}
-    template <unsigned n> explicit ConstData(const VectorGeneric<n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n> explicit ConstData(const VectorGeneric<n> &p): pointer(&p), size(n), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n,unsigned m> explicit ConstData(const TensorGeneric<n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n,unsigned m> explicit ConstData(const TensorGeneric<n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T,unsigned n> explicit ConstData(const VectorTyped<T,n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T,unsigned n> explicit ConstData(const VectorTyped<T,n> &p): pointer(&p), size(n), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T, unsigned n,unsigned m> explicit ConstData(const TensorTyped<T,n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T, unsigned n,unsigned m> explicit ConstData(const TensorTyped<T,n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(T)), type(getMPIType<T>()) {}
     template <typename T> explicit ConstData(const std::vector<T>&v) {
       ConstData d(v.data(),v.size());
       pointer=d.pointer;

src/tools/LoopUnroller.h

@@ -47,114 +47,89 @@ Implementation is made using template metaprogramming, that is:
 Here xxx is any of the methods of the class.
 
 */
-template<unsigned n>
+//the typename is the second parameter argument so that it can be deduced
+template<typename T, unsigned n>
 class LoopUnroller {
 public:
 /// Set to zero.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]=0.0;`
-  static void _zero(double*d);
+  constexpr static void _zero(T*d);
 /// Add v to d.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]+=v[i];`
-  static void _add(double*d,const double*v);
+  constexpr static void _add(T*d,const T*v);
 /// Subtract v from d.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]-=v[i];`
-  static void _sub(double*d,const double*v);
+  constexpr static void _sub(T*d,const T*v);
 /// Multiply d by s.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]*=s;`
-  static void _mul(double*d,const double s);
+  constexpr static void _mul(T*d,const T s);
 /// Set d to -v.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]=-v[i];`
-  static void _neg(double*d,const double*v);
+  constexpr static void _neg(T*d,const T*v);
 /// Squared modulo of d;
 /// Same as `r=0.0; for(unsigned i=0;i<n;i++) r+=d[i]*d[i]; return r;`
-  static double _sum2(const double*d);
+  constexpr static T _sum2(const T*d);
 /// Dot product of d and v
 /// Same as `r=0.0; for(unsigned i=0;i<n;i++) r+=d[i]*v[i]; return r;`
-  static double _dot(const double*d,const double*v);
+  constexpr static T _dot(const T*d,const T*v);
 };
 
-template<unsigned n>
-void LoopUnroller<n>::_zero(double*d) {
-  LoopUnroller<n-1>::_zero(d);
-  d[n-1]=0.0;
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_zero(T*d) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_zero(d);
+  }
+  d[n-1]=T(0);
 }
 
-template<>
-inline
-void LoopUnroller<1>::_zero(double*d) {
-  d[0]=0.0;
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_add(double*d,const double*a) {
-  LoopUnroller<n-1>::_add(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_add(T*d,const T*a) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_add(d,a);
+  }
   d[n-1]+=a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_add(double*d,const double*a) {
-  d[0]+=a[0];
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_sub(double*d,const double*a) {
-  LoopUnroller<n-1>::_sub(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_sub(T*d,const T*a) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_sub(d,a);
+  }
   d[n-1]-=a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_sub(double*d,const double*a) {
-  d[0]-=a[0];
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_mul(double*d,const double s) {
-  LoopUnroller<n-1>::_mul(d,s);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_mul(T*d,const T s) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_mul(d,s);
+  }
   d[n-1]*=s;
 }
 
-template<>
-inline
-void LoopUnroller<1>::_mul(double*d,const double s) {
-  d[0]*=s;
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_neg(double*d,const double*a ) {
-  LoopUnroller<n-1>::_neg(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_neg(T*d,const T*a ) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_neg(d,a);
+  }
   d[n-1]=-a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_neg(double*d,const double*a) {
-  d[0]=-a[0];
+template<typename T, unsigned n>
+constexpr T LoopUnroller<T,n>::_sum2(const T*d) {
+  if constexpr (n>1) {
+    return LoopUnroller<T,n-1>::_sum2(d)+d[n-1]*d[n-1];
+  } else {
+    return d[0]*d[0];
+  }
 }
-
-template<unsigned n>
-double LoopUnroller<n>::_sum2(const double*d) {
-  return LoopUnroller<n-1>::_sum2(d)+d[n-1]*d[n-1];
-}
-
-template<>
-inline
-double LoopUnroller<1>::_sum2(const double*d) {
-  return d[0]*d[0];
-}
-
-template<unsigned n>
-double LoopUnroller<n>::_dot(const double*d,const double*v) {
-  return LoopUnroller<n-1>::_dot(d,v)+d[n-1]*v[n-1];
-}
-
-template<>
-inline
-double LoopUnroller<1>::_dot(const double*d,const double*v) {
-  return d[0]*v[0];
-}
-
+template<typename T, unsigned n>
+constexpr T LoopUnroller<T,n>::_dot(const T*d,const T*v) {
+  if constexpr (n>1) {
+    return LoopUnroller<T,n-1>::_dot(d,v)+d[n-1]*v[n-1];
+  } else {
+    return d[0]*v[0];
+  }
 }
+} //PLMD
 
-#endif
+#endif //__PLUMED_tools_LoopUnroller_h

src/tools/MatrixSquareBracketsAccess.h

@@ -73,60 +73,60 @@ class MatrixSquareBracketsAccess {
     // the user should not manipulate it directly
     const MatrixSquareBracketsAccess& t;
     const I i;
-    Const_row(const MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
+    constexpr Const_row(const MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
   public:
     /// access element
-    const C & operator[] (J j)const;
+    constexpr const C & operator[] (J j)const;
   };
 /// Small utility class which just contains a pointer to the T and the row number
   class Row {
     friend class MatrixSquareBracketsAccess; // this so as to allow only T to instantiate Const_row
     // the user should not manipulate it directly
     MatrixSquareBracketsAccess& t;
     const I i;
-    Row(MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
+    constexpr Row(MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
   public:
     /// access element
-    C & operator[] (J j);
+    constexpr C & operator[] (J j);
   };
 public:
 /// access element (with [][] syntax)
-  Row operator[] (I i);
+  constexpr Row operator[] (I i);
 /// access element (with [][] syntax)
-  Const_row operator[] (I i)const;
+  constexpr Const_row operator[] (I i)const;
 };
 
 template<class T,class C,class I,class J>
-MatrixSquareBracketsAccess<T,C,I,J>::Const_row::Const_row(const MatrixSquareBracketsAccess&t,I i):
+constexpr MatrixSquareBracketsAccess<T,C,I,J>::Const_row::Const_row(const MatrixSquareBracketsAccess&t,I i):
   t(t),i(i) {}
 
 template<class T,class C,class I,class J>
-MatrixSquareBracketsAccess<T,C,I,J>::Row::Row(MatrixSquareBracketsAccess&t,I i):
+constexpr MatrixSquareBracketsAccess<T,C,I,J>::Row::Row(MatrixSquareBracketsAccess&t,I i):
   t(t),i(i) {}
 
 template<class T,class C,class I,class J>
-const C & MatrixSquareBracketsAccess<T,C,I,J>::Const_row::operator[] (J j)const {
+constexpr const C & MatrixSquareBracketsAccess<T,C,I,J>::Const_row::operator[] (J j)const {
 // This appears as a reference to a temporary object
 // however, in reality we know it is a reference to an object that is stored in the
 // t object. We thus suppress the warning raised by cppcheck
   return (*static_cast<const T*>(&t))(i,j);
 }
 
 template<class T,class C,class I,class J>
-C & MatrixSquareBracketsAccess<T,C,I,J>::Row::operator[] (J j) {
+constexpr C & MatrixSquareBracketsAccess<T,C,I,J>::Row::operator[] (J j) {
 // This appears as a reference to a temporary object
 // however, in reality we know it is a reference to an object that is stored in the
 // t object. We thus suppress the warning raised by cppcheck
   return (*static_cast<T*>(&t))(i,j);
 }
 
 template<class T,class C,class I,class J>
-typename MatrixSquareBracketsAccess<T,C,I,J>::Row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i) {
+constexpr typename MatrixSquareBracketsAccess<T,C,I,J>::Row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i) {
   return Row(*this,i);
 }
 
 template<class T,class C,class I,class J>
-typename MatrixSquareBracketsAccess<T,C,I,J>::Const_row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i)const {
+constexpr typename MatrixSquareBracketsAccess<T,C,I,J>::Const_row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i)const {
   return Const_row(*this,i);
 }