From 5ffaf0a1bb0f0b47d57d0f24ee1134659775dacb Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 20 Jul 2018 22:57:39 -0400
Subject: added ising example to cpp collection

---
 lib/vector.h | 98 +++++++++++++++++++++++++++++++-----------------------------
 1 file changed, 51 insertions(+), 47 deletions(-)

(limited to 'lib/vector.h')

diff --git a/lib/vector.h b/lib/vector.h
index a1084e2..c478618 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -6,16 +6,48 @@
 
 #include "types.h"
 
+/* The following is the minimum definition of a spin class.
+ *
+ * The class must contain an M_t and an F_t for holding the sum of an
+ * integer number of spins and a double-weighted number of spins,
+ * respectively.
+ *
+ * void init(X_t *p);
+ * void free_spin(X_t p);
+ * X_t copy(X_t x);
+ * void add(M_t *x1, int factor, X_t x2);
+ * void add(F_t *x1, double factor, X_t x2);
+ * M_t scalar_multiple(int factor, X_t x);
+ * double norm_squared(F_t x);
+ * void write_magnetization(M_t M, FILE *outfile);
+ *
+ */
+
 template <q_t q, class T>
-struct vector_t { T *x; };
+class vector_t {
+  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;
@@ -29,46 +61,15 @@ vector_t <q, T> copy (vector_t <q, T> v) {
   return v_copy;
 }
 
-template <q_t q, class T>
-void add (vector_t <q, T> *v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void scalar_add (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] += a * v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void scalar_subtract (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] -= a * v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
-  for (q_t i = 0; i < q; i++) {
-    v1->x[i] -= v2.x[i];
-  }
-}
-
-template <q_t q, class T>
-T norm_squared (vector_t <q, T> v) {
-  T tmp = 0;
+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++) {
-    tmp += v.x[i] * v.x[i];
+    v1->x[i] += (U)(a * v2.x[i]);
   }
-
-  return tmp;
 }
 
 template <q_t q, class T>
-vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
+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++) {
@@ -79,19 +80,14 @@ vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
 }
 
 template <q_t q, class T>
-T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
-  T prod = 0;
+double norm_squared (vector_t <q, T> v) {
+  double tmp = 0;
 
   for (q_t i = 0; i < q; i++) {
-    prod += v1.x[i] * v2.x[i];
+    tmp += pow(v.x[i], 2);
   }
 
-  return prod;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
-  free(v.x);
+  return tmp;
 }
 
 template <q_t q, class T>
@@ -99,6 +95,8 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) {
   fwrite(M.x, sizeof(T), q, outfile);
 }
 
+// below functions and definitions are unnecessary for wolff.h but useful.
+
 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++) {
@@ -108,8 +106,14 @@ void write_magnetization(vector_t <q, double> M, FILE *outfile) {
 }
 
 template <q_t q, class T>
-double correlation_component(vector_t <q, T> v) {
-  return (double)v.x[0];
+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];
+  }
+
+  return prod;
 }
 
 template <q_t q, class T>
-- 
cgit v1.2.3-70-g09d2