partially class-ified, ising and On work but potts and height do not

1 files changed, 57 insertions, 66 deletions

diff --git a/lib/vector.h b/lib/vector.h
index 2fe6ab8..d5ced03 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -3,6 +3,7 @@
 
 #include <stdlib.h>
 #include <cmath>
+#include <array>
 
 #include "types.h"
 
@@ -24,78 +25,67 @@
  */
 
 template <q_t q, class T>
-class vector_t {
+class vector_t : public std::array<T, q> {
   public: 
-    T *x;
 
     // M_t needs to hold the sum of nv spins
     typedef vector_t <q, T> M_t;
 
     // F_t needs to hold the double-weighted sum of spins
     typedef vector_t <q, double> F_t;
-};
-
-template <q_t q, class T>
-void init(vector_t <q, T> *ptr) {
-  ptr->x = (T *)calloc(q, sizeof(T));
-
-  // initialize a unit vector
-  ptr->x[0] = (T)1;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
-  free(v.x);
-}
-
-template <q_t q, class T>
-vector_t <q, T> copy (vector_t <q, T> v) {
-  vector_t <q, T> v_copy;
- 
- v_copy.x = (T *)calloc(q, sizeof(T));
 
-  for (q_t i = 0; i < q; i++) {
-    v_copy.x[i] = v.x[i];
-  }
-
-  return v_copy;
-}
-
-template <q_t q, class T, class U, class V>
-void add(vector_t<q, U> *v1, V a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += (U)(a * v2.x[i]);
-  }
-}
-
-template <q_t q, class T>
-vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) {
-  vector_t <q, T> multiple;
-  multiple.x = (T *)malloc(q * sizeof(T));
-  for (q_t i = 0; i < q; i++) {
-    multiple.x[i] = a * v.x[i];
-  }
-
-  return multiple;
-}
-
-template <q_t q, class T>
-vector_t <q, T> scalar_multiple(double a, vector_t <q, T> v) {
-  vector_t <q, T> multiple;
-  multiple.x = (T *)malloc(q * sizeof(T));
-  for (q_t i = 0; i < q; i++) {
-    multiple.x[i] = a * v.x[i];
-  }
+    vector_t() {
+      this->fill((T)0);
+      (*this)[1] = (T)1;
+    }
+
+    vector_t(const T *x) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] = x[i];
+      }
+    }
+
+    template <class U>
+    inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] += (U)v[i];
+      }
+      return *this;
+    }
+
+    template <class U>
+    inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] -= (U)v[i];
+      }
+      return *this;
+    }
+
+    inline vector_t<q, T> operator*(v_t x) const {
+      vector_t<q, T> result;
+      for (q_t i = 0; i < q; i++) {
+        result[i] = x * (*this)[i];
+      }
+
+      return result;
+    }
+
+    inline vector_t<q, double> operator*(double x) const {
+      vector_t<q, double> result;
+      for (q_t i = 0; i < q; i++) {
+        result[i] = x * (*this)[i];
+      }
+
+      return result;
+    }
+};
 
-  return multiple;
-}
 
-template <q_t q, class T>
-double norm_squared (vector_t <q, T> v) {
+template<q_t q>
+double norm_squared(vector_t<q, double> v) {
   double tmp = 0;
-
-  for (q_t i = 0; i < q; i++) {
-    tmp += pow(v.x[i], 2);
+  for (double &x : v) {
+    tmp += pow(x, 2);
   }
 
   return tmp;
@@ -103,7 +93,9 @@ double norm_squared (vector_t <q, T> v) {
 
 template <q_t q, class T>
 void write_magnetization(vector_t <q, T> M, FILE *outfile) {
-  fwrite(M.x, sizeof(T), q, outfile);
+  for (q_t i = 0; i < q; i++) {
+    fwrite(&(M[i]), sizeof(T), q, outfile);
+  }
 }
 
 // below functions and definitions are unnecessary for wolff.h but useful.
@@ -111,7 +103,7 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) {
 template <q_t q> // save some space and don't write whole doubles
 void write_magnetization(vector_t <q, double> M, FILE *outfile) {
   for (q_t i = 0; i < q; i++) {
-    float M_tmp = (float)M.x[i];
+    float M_tmp = (float)M[i];
     fwrite(&M_tmp, sizeof(float), 1, outfile);
   }
 }
@@ -121,7 +113,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
   T prod = 0;
 
   for (q_t i = 0; i < q; i++) {
-    prod += v1.x[i] * v2.x[i];
+    prod += v1[i] * v2[i];
   }
 
   return prod;
@@ -129,8 +121,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
 
 template <q_t q, class T>
 double H_vector(vector_t <q, T> v1, T *H) {
-  vector_t <q, T> H_vec;
-  H_vec.x = H;
+  vector_t <q, T> H_vec(H);
   return (double)(dot <q, T> (v1, H_vec));
 }