From 6cd2ea8fbdce43583d86fc195f0a62a8882cbeb6 Mon Sep 17 00:00:00 2001
From: Benjamin Christoffersen <boennecd@gmail.com>
Date: Sat, 30 Mar 2019 17:23:03 +0100
Subject: [PATCH] write custom matrix copy function

---
 NEWS.md              |  2 +-
 src/parallel_glm.cpp | 16 ++++++++++++++--
 2 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 66701fc..86416c6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,5 @@
 # pargm 0.1.3
-* Fix bug found by Valgrind.
+* Fix bug found with Valgrind.
 
 # pargm 0.1.2
 * Minor changes in implementation.
diff --git a/src/parallel_glm.cpp b/src/parallel_glm.cpp
index eb42dc5..3841d3f 100644
--- a/src/parallel_glm.cpp
+++ b/src/parallel_glm.cpp
@@ -5,7 +5,6 @@
 #include "R_BLAS_LAPACK.h"
 #include <memory>
 
-/* #define PARGLM_PROF */
 #ifdef PARGLM_PROF
 #include <gperftools/profiler.h>
 #include <iostream>
@@ -66,6 +65,19 @@ const double double_one = 1., double_zero = 0.;
 const int int_one = 1;
 char char_N = 'N', char_U = 'U', char_T = 'T';
 
+inline void inplace_copy
+  (arma::mat &X, const arma::mat &Y, const arma::uword start,
+   const arma::uword end)
+  {
+    double *x = X.begin();
+    const double *y = Y.begin() + start;
+
+    std::size_t n_ele = end - start + 1L;
+    for(unsigned int i = 0; i < X.n_cols;
+        ++i, x += X.n_rows, y += Y.n_rows)
+      std::memcpy(x, y, n_ele * sizeof(double));
+  }
+
 /* Class to fit glm using QR updated in chunks */
 class parallelglm_class_QR {
   using uword = arma::uword;
@@ -126,7 +138,7 @@ class parallelglm_class_QR {
 
       const arma::uword p = data.X.n_cols;
       arma::mat X(data.X_work_mem.get() + i_start * p, n, p, false, true);
-      X = data.X.rows(i_start, i_end);
+      inplace_copy(X, data.X, i_start, i_end);
       X.each_col() %= w;
       z %= w;