Dynare++: rollback the std::shared_ptr change in Vector and ConstVector

The performance cost is too high to warrant the change.

dynare++/kord/tests.cc

@@ -17,8 +17,7 @@ struct Rand
 ConstTwoDMatrix
 make_matrix(int rows, int cols, const double *p)
 {
-  return ConstTwoDMatrix{rows, cols,
+  return ConstTwoDMatrix{rows, cols, ConstVector{p, rows*cols}};
-      ConstVector{std::shared_ptr<const double>{p, [](const double *arr) {}}, rows*cols}};
 }
 
 void

dynare++/sylv/cc/Vector.cc

@@ -13,15 +13,15 @@
 #include <iomanip>
 
 Vector::Vector(const Vector &v)
-  : len(v.len), data{new double[len], [](double *arr) { delete[] arr; }}
+  : len(v.len), data{new double[len]}
 {
-  copy(v.base(), v.s);
+  copy(v.data, v.s);
 }
 
 Vector::Vector(const ConstVector &v)
-  : len(v.len), data{new double[len], [](double *arr) { delete[] arr; }}
+  : len(v.len), data{new double[len]}
 {
-  copy(v.base(), v.s);
+  copy(v.data, v.s);
 }
 
 Vector &
@@ -34,12 +34,12 @@ Vector::operator=(const Vector &v)
     throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
 
   if (s == v.s
-      && (base() <= v.base() && v.base() < base()+len*s
+      && (data <= v.data && v.data < data+len*s
-          || v.base() <= base() && base() < v.base()+v.len*v.s)
+          || v.data <= data && data < v.data+v.len*v.s)
-      && (base()-v.base()) % s == 0)
+      && (data-v.data) % s == 0)
     throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
 
-  copy(v.base(), v.s);
+  copy(v.data, v.s);
   return *this;
 }
 
@@ -48,7 +48,7 @@ Vector::operator=(Vector &&v)
 {
   if (v.len != len)
     throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
-  copy(v.base(), v.s);
+  copy(v.data, v.s);
   return *this;
 }
 
@@ -58,12 +58,12 @@ Vector::operator=(const ConstVector &v)
   if (v.len != len)
     throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
   if (s == v.s
-      && (base() <= v.base() && v.base() < base()+len*s
+      && (data <= v.data && v.data < data+len*s
-          || v.base() <= base() && base() < v.base()+v.len*v.s)
+          || v.data <= data && data < v.data+v.len*v.s)
-      && (base()-v.base()) % s == 0)
+      && (data-v.data) % s == 0)
     throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
 
-  copy(v.base(), v.s);
+  copy(v.data, v.s);
   return *this;
 }
 
@@ -73,39 +73,39 @@ Vector::copy(const double *d, int inc)
   blas_int n = len;
   blas_int incy = s;
   blas_int inc2 = inc;
-  dcopy(&n, d, &inc2, base(), &incy);
+  dcopy(&n, d, &inc2, data, &incy);
 }
 
 Vector::Vector(Vector &v, int off_arg, int l)
-  : len(l), off{v.off+off_arg*v.s}, s(v.s), data{v.data}
+  : len(l), s(v.s), data{v.data+off_arg*v.s}, destroy{false}
 {
   if (off_arg < 0 || off_arg + len > v.len)
     throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
 }
 
 Vector::Vector(const Vector &v, int off_arg, int l)
-  : len(l), data{new double[len], [](double *arr) { delete[] arr; }}
+  : len(l), data{new double[len]}
 {
   if (off_arg < 0 || off_arg + len > v.len)
     throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
-  copy(v.base()+off_arg*v.s, v.s);
+  copy(v.data+off_arg*v.s, v.s);
 }
 
 Vector::Vector(Vector &v, int off_arg, int skip, int l)
-  : len(l), off{v.off+off_arg*v.s}, s(v.s*skip), data{v.data}
+  : len(l), s(v.s*skip), data{v.data+off_arg*v.s}, destroy{false}
 {
 }
 
 Vector::Vector(const Vector &v, int off_arg, int skip, int l)
-  : len(l), data{new double[len], [](double *arr) { delete[] arr; }}
+  : len(l), data{new double[len]}
 {
-  copy(v.base()+off_arg*v.s, v.s*skip);
+  copy(v.data+off_arg*v.s, v.s*skip);
 }
 
 #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 Vector::Vector(mxArray *p)
   : len{static_cast<int>(mxGetNumberOfElements(p))},
-    data{std::shared_ptr<double>{mxGetPr(p), [](double *arr) { }}}
+    data{mxGetPr(p)}, destroy{false}
 {
   if (!mxIsDouble(p))
     throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles.");
@@ -152,7 +152,7 @@ void
 Vector::zeros()
 {
   if (s == 1)
-    std::fill_n(base(), len, 0.0);
+    std::fill_n(data, len, 0.0);
   else
     for (int i = 0; i < len; i++)
       operator[](i) = 0.0;
@@ -192,7 +192,7 @@ Vector::add(double r, const ConstVector &v)
   blas_int n = len;
   blas_int incx = v.s;
   blas_int incy = s;
-  daxpy(&n, &r, v.base(), &incx, base(), &incy);
+  daxpy(&n, &r, v.data, &incx, data, &incy);
 }
 
 void
@@ -207,7 +207,7 @@ Vector::addComplex(const std::complex<double> &z, const ConstVector &v)
   blas_int n = len/2;
   blas_int incx = v.s;
   blas_int incy = s;
-  zaxpy(&n, reinterpret_cast<const double(&)[2]>(z), v.base(), &incx, base(), &incy);
+  zaxpy(&n, reinterpret_cast<const double(&)[2]>(z), v.data, &incx, data, &incy);
 }
 
 void
@@ -215,7 +215,7 @@ Vector::mult(double r)
 {
   blas_int n = len;
   blas_int incx = s;
-  dscal(&n, &r, base(), &incx);
+  dscal(&n, &r, data, &incx);
 }
 
 void
@@ -283,31 +283,31 @@ Vector::print() const
 }
 
 ConstVector::ConstVector(const Vector &v)
-  : len{v.len}, off{v.off}, s{v.s}, data{v.data}
+  : len{v.len}, s{v.s}, data{v.data}
 {
 }
 
 ConstVector::ConstVector(const ConstVector &v, int off_arg, int l)
-  : len{l}, off{v.off+off_arg*v.s}, s{v.s}, data{v.data}
+  : len{l}, s{v.s}, data{v.data+off_arg*v.s}
 {
   if (off_arg < 0 || off_arg + len > v.len)
     throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
 }
 
 ConstVector::ConstVector(const ConstVector &v, int off_arg, int skip, int l)
-  : len(l), off{v.off+off_arg*v.s}, s{v.s*skip}, data{v.data}
+  : len(l), s{v.s*skip}, data{v.data+off_arg*v.s}
 {
 }
 
-ConstVector::ConstVector(std::shared_ptr<const double> d, int skip, int l)
+ConstVector::ConstVector(const double *d, int skip, int l)
-  : len{l}, s{skip}, data{std::move(d)}
+  : len{l}, s{skip}, data{d}
 {
 }
 
 #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 ConstVector::ConstVector(const mxArray *p)
   : len{static_cast<int>(mxGetNumberOfElements(p))},
-    data{std::shared_ptr<const double>{mxGetPr(p), [](const double *arr) { }}}
+    data{mxGetPr(p)}
 {
   if (!mxIsDouble(p))
     throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles.");
@@ -376,7 +376,7 @@ ConstVector::dot(const ConstVector &y) const
   blas_int n = len;
   blas_int incx = s;
   blas_int incy = y.s;
-  return ddot(&n, base(), &incx, y.base(), &incy);
+  return ddot(&n, data, &incx, y.data, &incy);
 }
 
 bool

dynare++/sylv/cc/Vector.hh

@@ -9,7 +9,6 @@
  * to avoid running virtual method invokation mechanism. Some
  * members, and methods are thus duplicated */
 
-#include <memory>
 #include <complex>
 
 #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
@@ -24,21 +23,30 @@ class Vector
   friend class ConstVector;
 protected:
   int len{0};
-  int off{0}; // offset to double* pointer
   int s{1}; // stride (also called "skip" in some places)
-  std::shared_ptr<double> data;
+  double *data;
+  bool destroy{true};
 public:
-  Vector() = default;
+  Vector() : data{nullptr}, destroy{false}
-  Vector(int l) : len{l}, data{new double[l], [](double *arr) { delete[] arr; }}
+  {
+  }
+  Vector(int l) : len{l}, data{new double[l]}
+  {
+  }
+  Vector(Vector &v) : len{v.len}, s{v.s}, data{v.data}, destroy{false}
   {
   }
-  Vector(Vector &v) = default;
   Vector(const Vector &v);
-  Vector(Vector &&v) = default;
+  Vector(Vector &&v) : len{v.len}, s{v.s}, data{v.data}, destroy{v.destroy}
+  {
+    v.len = 0;
+    v.data = nullptr;
+    v.destroy = false;
+  }
   // We don't want implict conversion from ConstVector, since it's expensive
   explicit Vector(const ConstVector &v);
-  Vector(std::shared_ptr<double> d, int l)
+  Vector(double *d, int l)
-    : len(l), data{std::move(d)}
+    : len(l), data{d}, destroy{false}
   {
   }
   Vector(Vector &v, int off_arg, int l);
@@ -54,22 +62,22 @@ public:
   double &
   operator[](int i)
   {
-    return data.get()[off+s*i];
+    return data[s*i];
   }
   const double &
   operator[](int i) const
   {
-    return data.get()[off+s*i];
+    return data[s*i];
   }
   const double *
   base() const
   {
-    return data.get() + off;
+    return data;
   }
   double *
   base()
   {
-    return data.get() + off;
+    return data;
   }
   int
   length() const
@@ -91,7 +99,11 @@ public:
   bool operator>(const Vector &y) const;
   bool operator>=(const Vector &y) const;
 
-  virtual ~Vector() = default;
+  virtual ~Vector()
+  {
+    if (destroy)
+      delete[] data;
+  }
   void zeros();
   void nans();
   void infs();
@@ -142,20 +154,19 @@ class ConstVector
   friend class Vector;
 protected:
   int len;
-  int off{0}; // offset to double* pointer
   int s{1}; // stride (also called "skip" in some places)
-  std::shared_ptr<const double> data;
+  const double *data;
 public:
   // Implicit conversion from Vector is ok, since it's cheap
   ConstVector(const Vector &v);
   ConstVector(const ConstVector &v) = default;
   ConstVector(ConstVector &&v) = default;
-  ConstVector(std::shared_ptr<const double> d, int l) : len{l}, data{std::move(d)}
+  ConstVector(const double *d, int l) : len{l}, data{d}
   {
   }
   ConstVector(const ConstVector &v, int off_arg, int l);
   ConstVector(const ConstVector &v, int off_arg, int skip, int l);
-  ConstVector(std::shared_ptr<const double> d, int skip, int l);
+  ConstVector(const double *d, int skip, int l);
 #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit ConstVector(const mxArray *p);
 #endif
@@ -165,12 +176,12 @@ public:
   const double &
   operator[](int i) const
   {
-    return data.get()[off+s*i];
+    return data[s*i];
   }
   const double *
   base() const
   {
-    return data.get() + off;
+    return data;
   }
   int
   length() const

dynare++/sylv/testing/MMMatrix.cc

@@ -22,12 +22,12 @@ MMMatrixIn::MMMatrixIn(const std::string &fname)
   if (fd.fail())
     throw MMException("Couldn't parse rows and cols\n");
   // read in data
-  data = std::shared_ptr<double>(static_cast<double *>(operator new[](rows*cols*sizeof(double))), [](double *arr) { operator delete[](static_cast<void *>(arr)); });
   int len = rows*cols;
+  data.resize(len);
   int i = 0;
   while (!fd.eof() && i < len)
     {
-      fd >> data.get()[i];
+      fd >> data[i];
       if (fd.fail())
         throw MMException("Couldn't parse float number\n");
       i++;

dynare++/sylv/testing/MMMatrix.hh

@@ -9,7 +9,7 @@
 
 #include <string>
 #include <utility>
-#include <memory>
+#include <vector>
 
 class MMException
 {
@@ -27,16 +27,16 @@ public:
 
 class MMMatrixIn
 {
-  std::shared_ptr<double> data;
+  std::vector<double> data;
   int rows;
   int cols;
 public:
   MMMatrixIn(const std::string &fname);
   ~MMMatrixIn() = default;
   Vector
-  getData() const
+  getData()
   {
-    return Vector{data, size()};
+    return Vector{data.data(), size()};
   }
   int
   size() const