Typed vector and tensor (plumed#1092)

* templatized runtime fastpow

* templatized the type of LoopUnroller

* templatized the type of Vector

* templatized the type of Tensor

* constexpr everithing!

* solving some cpp checks

* rounded a few compatibilities between the new Vector and View

* Updating Vector and Tensor with conversion and types interactions

* added and checked the possibility of using ssyevr

* lowestEigenpairSym typing

* Changing the precision in rt-make-4

* adressed codecheck uninitialized return variable

---------

Co-authored-by: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>

Author: Iximiel

regtest/basic/rt-lowest-eigenpair/eigenValues_float.reference

@@ -16,8 +16,9 @@ using namespace PLMD;
 
 namespace {
 
-static Tensor4d makeMatrix(const Tensor & rr01) {
-  Tensor4d m;
+template<typename precision>
+static auto makeMatrix(const TensorT<precision> & rr01) {
+  TensorTyped<precision,4,4> m;
   m[0][0]= +rr01[0][0]+rr01[1][1]+rr01[2][2];
   m[1][1]= +rr01[0][0]-rr01[1][1]-rr01[2][2];
   m[2][2]= -rr01[0][0]+rr01[1][1]-rr01[2][2];
@@ -51,9 +52,9 @@ inline void DoNotOptimize(T& value) {
 }
 
 
-
-PLMD::Tensor randomRMSDmat(rng & gen, const double referencedistance=1.0) {
-  PLMD::Tensor toret;
+template <typename precision>
+auto randomRMSDmat(rng & gen, const double referencedistance=1.0) {
+  TensorT<precision> toret;
   constexpr size_t n=30;
 
   std::vector<double>  align(n);
@@ -86,16 +87,17 @@ PLMD::Tensor randomRMSDmat(rng & gen, const double referencedistance=1.0) {
 // second expensive loop: compute second moments wrt centers
   for(unsigned iat=0; iat<n; iat++) {
     double w=align[iat];
-    toret+=Tensor(positions[iat]-cpositions,reference[iat])*w;
+    PLMD::LoopUnroller<9>::_add(toret.data(),
+                                (Tensor(positions[iat]-cpositions,reference[iat])*w).data());
   }
   return toret;
 }
 
-
-int main() {
+template <typename T>
+void doTest(std::string suffix="") {
   std::cerr <<std::setprecision(9) << std::fixed;
 
-  std::vector<PLMD::Tensor> tests= {
+  std::vector<PLMD::TensorT<T>> tests= {
     // {1,1,1,1,1,1,1,1,1}, // 'wrong' eigenvector selected (due to the degeneracy)
     // {1,2,3,4,5,6,7,8,9}, // 'wrong' eigenvector selected (due to the degeneracy)
     {0,1,1,1,1,1,1,1,1},
@@ -125,14 +127,14 @@ int main() {
   gen.setSeed(123456789);
   double error=100.0;
   for(int i=0; i<400; i++) {
-    tests.push_back(randomRMSDmat(gen,error));
+    tests.push_back(randomRMSDmat<T>(gen,error));
     if(i%50) {
       error/=10.0;
     }
   }
-  std::ofstream matOut ("matrices");
-  std::ofstream out ("eigenValues");
-  std::ofstream outv ("eigenVectors");
+  std::ofstream matOut ("matrices"+suffix);
+  std::ofstream out ("eigenValues"+suffix);
+  std::ofstream outv ("eigenVectors"+suffix);
 // setting up an header for clarity
 
   matOut<<"# id:";
@@ -162,12 +164,12 @@ int main() {
     const auto& mat = tests[index];
     auto m=makeMatrix(mat);
 
-    Vector1d eigenval_ref;
-    TensorGeneric<1,4> eigenvec_ref;
+    VectorTyped<T,1> eigenval_ref;
+    TensorTyped<T,1,4> eigenvec_ref;
     diagMatSym(m, eigenval_ref, eigenvec_ref );
 
-    double eigenval;
-    Vector4d eigenvec;
+    T eigenval;
+    PLMD::VectorTyped<T,4> eigenvec;
     eigenval=lowestEigenpairSym(m,eigenvec);
 
     matOut <<std::setw(4)<<index <<":";
@@ -178,7 +180,12 @@ int main() {
     matOut<<"\n";
 
     out <<std::setw(4)<<index <<":";
-    out<<" "<<std::setw(10) <<eigenval_ref[0] - eigenval;
+    if constexpr (std::is_same_v<T,double>) {
+      out<<" "<<std::setw(10) <<eigenval_ref[0] - eigenval;
+    } else {
+      auto tmp = eigenval_ref[0] - eigenval;
+      out<<" "<<std::setw(10) <<((std::fabs(tmp) < 1e-4) ? T(0.0): tmp);
+    }
     out<<"\n";
 
     double modulo2_ref=0.0;
@@ -198,7 +205,7 @@ int main() {
 
   //Timings
   {
-    std::vector<PLMD::Tensor4d> matrices(tests.size());
+    std::vector<PLMD::TensorTyped<T,4,4>> matrices(tests.size());
     for(size_t i=0; i< tests.size(); ++i) {
       matrices[i] = makeMatrix(tests[i]);
     }
@@ -210,8 +217,8 @@ int main() {
         for (const auto & mat : matrices) {
           auto sww = sw.startStop("diagMatSym");
 
-          Vector1d eigenvals;
-          TensorGeneric<1,4> eigenvecs;
+          VectorTyped<T,1> eigenvals;
+          TensorTyped<T,1,4> eigenvecs;
 
           diagMatSym(mat, eigenvals, eigenvecs );
           DoNotOptimize(eigenvals);
@@ -224,16 +231,21 @@ int main() {
         for (const auto & mat : matrices) {
           auto sww = sw.startStop("compute_quaternion_from_K");
 
-          double eigenvals;
-          Vector4d eigenvecs;
+          VectorTyped<T,4> eigenvecs;
 
-          eigenvals=lowestEigenpairSym(mat, eigenvecs);
+          auto eigenvals=lowestEigenpairSym(mat, eigenvecs);
           DoNotOptimize(eigenvals);
           DoNotOptimize(eigenvecs);
         }
       }
     }
     std::cerr<<sw;
   }
+}
+int main () {
+  doTest<double>();
+  std::cerr << "\n\nFloating point\n\n";
+  doTest<float>("_float");
+
   return 0;
 }

regtest/basic/rt-lowest-eigenpair/eigenVectors_float.reference

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

regtest/basic/rt-lowest-eigenpair/main.cpp

@@ -2,42 +2,60 @@
 #include "plumed/tools/Matrix.h"
 #include "plumed/tools/Tensor.h"
 
-int main () {
+template <typename precision>
+void matrixTest (PLMD::OFile& out, const std::string& fmt) {
 
   // Define symmetric matrix
-  PLMD::Matrix<double> mat1(3,3); PLMD::OFile out; out.open("output");
-  mat1(0,0)=1.0; mat1(0,1)=0.2; mat1(0,2)=0.3;
-  mat1(1,0)=0.2; mat1(1,1)=0.2; mat1(1,2)=0.6;
-  mat1(2,0)=0.3; mat1(2,1)=0.6; mat1(2,2)=0.4;
+  PLMD::Matrix<precision> mat1(3,3);
+  mat1(0,0)=1.0;
+  mat1(0,1)=0.2;
+  mat1(0,2)=0.3;
+  mat1(1,0)=0.2;
+  mat1(1,1)=0.2;
+  mat1(1,2)=0.6;
+  mat1(2,0)=0.3;
+  mat1(2,1)=0.6;
+  mat1(2,2)=0.4;
 
   // Test diagonalize
-  std::vector<double> eigval(3); PLMD::Matrix<double> eigvec(3,3); 
-  diagMat( mat1, eigval, eigvec ); 
-  out<<"Eigenvalues "<<eigval[0]<<" "<<eigval[1]<<" "<<eigval[2]<<"\n";
+  std::vector<precision> eigval(3);
+  PLMD::Matrix<precision> eigvec(3,3);
+  diagMat( mat1, eigval, eigvec );
+  out.printf(("Eigenvalues "+fmt+" "+fmt+" "+fmt+"\n").c_str(),
+             eigval[0],eigval[1],eigval[2]);
   out<<"Eigenvectors : \n";
-  for(unsigned i=0;i<3;++i){
-      out<<eigvec(i,0)<<" "<<eigvec(i,1)<<" "<<eigvec(i,2)<<"\n";
+  for(unsigned i=0; i<3; ++i) {
+    out.printf((fmt+" "+fmt+" "+fmt+"\n").c_str(),
+               eigvec(i,0), eigvec(i,1), eigvec(i,2));
   }
 
   // Test inverse
   out<<"Inverse : \n";
-  PLMD::Matrix<double> inverse(3,3); Invert( mat1, inverse );
-  for(unsigned i=0;i<3;++i){ 
-      for(unsigned j=0;j<3;++j) out<<inverse(i,j)<<" ";
-      out<<"\n";
+  PLMD::Matrix<precision> inverse(3,3);
+  Invert( mat1, inverse );
+  for(unsigned i=0; i<3; ++i) {
+    for(unsigned j=0; j<3; ++j) {
+      out.printf((" "+fmt).c_str(),inverse(i,j));
+    }
+    out<<"\n";
   }
 
-  // Test pseudoinverse 
+  // Test pseudoinverse
   out<<"Pseudoinverse : \n";
-  PLMD::Matrix<double> mat(3,2);
-  mat(0,0)=0.1; mat(0,1)=0.2; 
-  mat(1,0)=0.3; mat(1,1)=0.5;
-  mat(2,0)=0.4; mat(2,1)=0.6;
-  PLMD::Matrix<double> pseu(2,3);
+  PLMD::Matrix<precision> mat(3,2);
+  mat(0,0)=0.1;
+  mat(0,1)=0.2;
+  mat(1,0)=0.3;
+  mat(1,1)=0.5;
+  mat(2,0)=0.4;
+  mat(2,1)=0.6;
+  PLMD::Matrix<precision> pseu(2,3);
   pseudoInvert( mat, pseu );
-  for(unsigned i=0;i<pseu.nrows();++i){
-     for(unsigned j=0;j<pseu.ncols();++j) out<<" "<<pseu(i,j);
-     out<<"\n";
+  for(unsigned i=0; i<pseu.nrows(); ++i) {
+    for(unsigned j=0; j<pseu.ncols(); ++j) {
+      out.printf((" "+fmt).c_str(),pseu(i,j));
+    }
+    out<<"\n";
   }
 
 
@@ -47,11 +65,13 @@ int main () {
 /// This would trigger an error before fix
 /// f1da0a9b3a13f904bd97d6f92a2fb5e1b6479ac0
   {
-    PLMD::TensorGeneric<4,4> mat;
-    PLMD::TensorGeneric<1,4> evec;
-    PLMD::VectorGeneric<4> eval_underlying;
+    PLMD::TensorTyped<precision,4,4> mat;
+    PLMD::TensorTyped<precision,1,4> evec;
+    PLMD::VectorTyped<precision,4> eval_underlying;
 
-    auto eval = new(&eval_underlying[0]) PLMD::VectorGeneric<1>;
+    //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::VectorTyped<precision,1>;
 
     mat[1][0]=mat[0][1]=3.0;
     mat[1][1]=5.0;
@@ -61,16 +81,16 @@ int main () {
 
     out<<"Test diagmat (identical eigenvalues)\n";
 /// Also eval_underlying[1] will be modified
-    out<<eval_underlying[0]<<" "<<eval_underlying[1];
-    out<<"\n";
+    out.printf((fmt+" "+fmt+"\n").c_str(),
+               eval_underlying[0], eval_underlying[1]);
   }
 
   {
-    PLMD::TensorGeneric<4,4> mat;
-    PLMD::TensorGeneric<1,4> evec;
-    PLMD::VectorGeneric<4> eval_underlying;
+    PLMD::TensorTyped<precision,4,4> mat;
+    PLMD::TensorTyped<precision,1,4> evec;
+    PLMD::VectorTyped<precision,4> eval_underlying;
 
-    auto eval = new(&eval_underlying[0]) PLMD::VectorGeneric<1>;
+    auto eval = new(&eval_underlying[0]) PLMD::VectorTyped<precision,1>;
 
     mat[1][0]=mat[0][1]=3.0;
     mat[1][1]=5.0;
@@ -79,18 +99,18 @@ int main () {
     PLMD::diagMatSym(mat,*eval,evec);
 
     out<<"Test diagmat (different eigenvalues)\n";
-    out<<eval_underlying[0]<<" "<<eval_underlying[1];
-    out<<"\n";
+    out.printf((fmt+" "+fmt+"\n").c_str(),
+               eval_underlying[0], eval_underlying[1]);
   }
 
 /// The following is to test the fact that
 /// dsyevr actually modifies the diagonalized matrix
 /// This would trigger an error before fix
 /// f1da0a9b3a13f904bd97d6f92a2fb5e1b6479ac0
   {
-    PLMD::TensorGeneric<4,4> mat;
-    PLMD::TensorGeneric<4,4> evec;
-    PLMD::VectorGeneric<4> eval;
+    PLMD::TensorTyped<precision,4,4> mat;
+    PLMD::TensorTyped<precision,4,4> evec;
+    PLMD::VectorTyped<precision,4> eval;
 
     mat[1][0]=mat[0][1]=3.0;
     mat[1][1]=5.0;
@@ -100,19 +120,32 @@ int main () {
 
     out<<"Test diagmat (is matrix modified)\n";
     out<<"Before:\n";
-    for(unsigned i=0;i<4;i++) {
-      for(unsigned j=0;j<4;j++) out<<" "<<mat[i][j];
+    for(unsigned i=0; i<4; i++) {
+      for(unsigned j=0; j<4; j++) {
+        out.printf((" "+fmt).c_str(),mat[i][j]);
+      }
       out<<"\n";
     }
 
     PLMD::diagMatSym(mat,eval,evec);
 
     out<<"After:\n";
-    for(unsigned i=0;i<4;i++) {
-      for(unsigned j=0;j<4;j++) out<<" "<<mat[i][j];
+    for(unsigned i=0; i<4; i++) {
+      for(unsigned j=0; j<4; j++) {
+        out.printf((" "+fmt).c_str(),mat[i][j]);
+      }
       out<<"\n";
     }
   }
+}
+int main () {
+  PLMD::OFile out;
+  out.open("output");
+  matrixTest<double>(out,"%8.6f");
   out.close();
+  PLMD::OFile out_float;
+  out_float.open("output_float");
+  matrixTest<float>(out_float,"%8.4f");
+  out_float.close();
   return 0;
 }

regtest/basic/rt-make-0/config

@@ -1,27 +1,27 @@
-Eigenvalues -0.309693 0.593857 1.31584
+Eigenvalues -0.309693 0.593857 1.315836
 Eigenvectors : 
-0.0342772 0.756017 -0.653654
+0.034277 0.756017 -0.653654
 0.659404 -0.508584 -0.553651
 0.751007 0.412044 0.515953
 Inverse : 
-1.15702 -0.413223 -0.247934 
--0.413223 -1.28099 2.2314 
--0.247934 2.2314 -0.661157 
+ 1.157025 -0.413223 -0.247934
+ -0.413223 -1.280992 2.231405
+ -0.247934 2.231405 -0.661157
 Pseudoinverse : 
- -18.8889 -11.1111 15.5556
- 12.2222 7.77778 -8.88889
+ -18.888889 -11.111111 15.555556
+ 12.222222 7.777778 -8.888889
 Test diagmat (identical eigenvalues)
--1.40512 0
+-1.405125 0.000000
 Test diagmat (different eigenvalues)
--1.40512 0
+-1.405125 0.000000
 Test diagmat (is matrix modified)
 Before:
- 0 3 0 -1
- 3 5 0 0
- 0 0 0 3
- -1 0 3 6
+ 0.000000 3.000000 0.000000 -1.000000
+ 3.000000 5.000000 0.000000 0.000000
+ 0.000000 0.000000 0.000000 3.000000
+ -1.000000 0.000000 3.000000 6.000000
 After:
- 0 3 0 -1
- 3 5 0 0
- 0 0 0 3
- -1 0 3 6
+ 0.000000 3.000000 0.000000 -1.000000
+ 3.000000 5.000000 0.000000 0.000000
+ 0.000000 0.000000 0.000000 3.000000
+ -1.000000 0.000000 3.000000 6.000000