From 3039e2c56532d4de852848df159df3041ee4bb4b Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 3 Feb 2020 16:07:45 -0500
Subject: Givens rotations now use single vector for storage instead of
 allocating one each construction

---
 hadamard_mcmc.hpp | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index 218f8d1..e64eacf 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -53,18 +53,16 @@ private:
   unsigned axis_2;
   double Δθ;
 
-  std::vector<double> rows;
-
 public:
   Givens(Orthogonal& m, bool t, unsigned a1, unsigned a2, double θ0, Rng& rng)
-      : m(m), rows(2 * m.size()) {
+      : m(m) {
     transpose = t;
     axis_1 = a1;
     axis_2 = a2;
     Δθ = rng.uniform(-θ0, θ0);
   }
 
-  Givens(Orthogonal& m, double θ0, Rng& rng) : m(m), rows(2 * m.size()) {
+  Givens(Orthogonal& m, double θ0, Rng& rng) : m(m) {
     Δθ = rng.uniform(-θ0, θ0);
     unsigned axis1axis2 = rng.uniform((unsigned)0, m.size() * (m.size() - 1) - 1);
 
@@ -76,7 +74,7 @@ public:
     }
   }
 
-  double tryRotation() {
+  double tryRotation(std::vector<double>& rows) const {
     double ΔE = 0;
     double c = cos(Δθ);
     double s = sin(Δθ);
@@ -106,7 +104,7 @@ public:
     return ΔE;
   }
 
-  void acceptRotation() const {
+  void acceptRotation(std::vector<double>& rows) const {
     for (unsigned i = 0; i < m.size(); i++) {
       if (transpose) {
         m(i, axis_1) = rows[i];
@@ -133,6 +131,7 @@ typedef enum {
 
 class MCMC {
 private:
+  std::vector<double> row_storage;
   Measurement& A;
   double θ0;
 
@@ -142,19 +141,19 @@ public:
   double E;
   Orthogonal M;
 
-  MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0) {
+  MCMC(unsigned n, double β0, Measurement& A) : A(A), M(n), β(β0), row_storage(2 * n) {
     θ0 = M_PI;
     E = M.energy();
   }
 
   bool step(Givens& g, bool dry = false) {
-    double ΔE = g.tryRotation();
+    double ΔE = g.tryRotation(row_storage);
 
     bool accepted = ΔE < 0 || exp(-β * ΔE) > rng.uniform((double)0.0, 1.0);
 
     if (accepted) {
       E += ΔE;
-      g.acceptRotation();
+      g.acceptRotation(row_storage);
     }
 
     if (!dry) {
-- 
cgit v1.2.3-54-g00ecf