diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc
index e81b7098213f9fd881cdd596adf42720b645b2c3..0bc2c70ff76f669356415a2e64fb914cfec86f71 100644
--- a/dynare++/extern/matlab/dynare_simul.cc
+++ b/dynare++/extern/matlab/dynare_simul.cc
@@ -110,8 +110,8 @@ extern "C" {
           dr(pol, PartitionY(nstat, npred, nboth, nforw),
              nexog, ConstVector{ysteady});
         // form the shock realization
-        TwoDMatrix shocks_mat(nexog, nper, ConstVector{shocks});
-        TwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov});
+        ConstTwoDMatrix shocks_mat(nexog, nper, ConstVector{shocks});
+        ConstTwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov});
         GenShockRealization sr(vcov_mat, shocks_mat, seed);
         // simulate and copy the results
         TwoDMatrix *res_mat
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index f61d937b5a21ca2e7eafc2947f055f3ac95ad337..aade85af6c59ca97c05928fe0ea81dc3c0c2cd9a 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -560,7 +560,7 @@ RTSimulationWorker::operator()()
    not work for semidefinite matrices. */
 
 void
-RandomShockRealization::choleskyFactor(const TwoDMatrix &v)
+RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
 {
   factor = v;
   lapack_int rows = factor.nrows(), lda = factor.getLD();
@@ -578,7 +578,7 @@ RandomShockRealization::choleskyFactor(const TwoDMatrix &v)
    decomposition. It works for semidifinite matrices. */
 
 void
-RandomShockRealization::schurFactor(const TwoDMatrix &v)
+RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
 {
   SymSchurDecomp ssd(v);
   ssd.getFactor(factor);
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 5b9fa80da32d83052feb36da55995d73cc6b6cbd..72f0081b0f24174dc9e8716e5758e93814a3e018 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -990,7 +990,7 @@ protected:
   MersenneTwister mtwister;
   TwoDMatrix factor;
 public:
-  RandomShockRealization(const TwoDMatrix &v, unsigned int iseed)
+  RandomShockRealization(const ConstTwoDMatrix &v, unsigned int iseed)
     : mtwister(iseed), factor(v.nrows(), v.nrows())
   {
     schurFactor(v);
@@ -1008,8 +1008,8 @@ public:
     return factor.nrows();
   }
 protected:
-  void choleskyFactor(const TwoDMatrix &v);
-  void schurFactor(const TwoDMatrix &v);
+  void choleskyFactor(const ConstTwoDMatrix &v);
+  void schurFactor(const ConstTwoDMatrix &v);
 };
 
 /* This is just a matrix of finite numbers. It can be constructed from
@@ -1019,10 +1019,6 @@ class ExplicitShockRealization : virtual public ShockRealization
 {
   TwoDMatrix shocks;
 public:
-  ExplicitShockRealization(const TwoDMatrix &sh)
-    : shocks(sh)
-  {
-  }
   ExplicitShockRealization(const ConstTwoDMatrix &sh)
     : shocks(sh)
   {
@@ -1063,7 +1059,7 @@ public:
 class GenShockRealization : public RandomShockRealization, public ExplicitShockRealization
 {
 public:
-  GenShockRealization(const TwoDMatrix &v, const TwoDMatrix &sh, int seed)
+  GenShockRealization(const ConstTwoDMatrix &v, const ConstTwoDMatrix &sh, int seed)
     : RandomShockRealization(v, seed), ExplicitShockRealization(sh)
   {
     KORD_RAISE_IF(sh.nrows() != v.nrows() || v.nrows() != v.ncols(),
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index ee2d4633ebed7e51d0c42e6b56839573d161c744..f721417ecded0bb608ae7373fa4c0d376a39f7fb 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -13,10 +13,10 @@ struct Rand
   static bool discrete(double prob); // answers true with given probability
 };
 
-TwoDMatrix
+ConstTwoDMatrix
 make_matrix(int rows, int cols, const double *p)
 {
-  return TwoDMatrix{rows, cols,
+  return ConstTwoDMatrix{rows, cols,
       ConstVector{std::shared_ptr<const double>{p, [](const double *arr) {}}, rows*cols}};
 }
 
diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc
index 7ba83d783ae6a0c0749918238a1a4da5cb66cbc3..20041cfe4d6880f98b5b4c1042a15dc7a433d78c 100644
--- a/dynare++/sylv/cc/GeneralMatrix.cc
+++ b/dynare++/sylv/cc/GeneralMatrix.cc
@@ -80,6 +80,16 @@ GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const char *dum1,
   gemm("T", a, "T", b, 1.0, 0.0);
 }
 
+GeneralMatrix &
+GeneralMatrix::operator=(const ConstGeneralMatrix &m)
+{
+  data = m.data;
+  rows = m.rows;
+  cols = m.cols;
+  ld = m.ld;
+  return *this;
+}
+
 Vector
 GeneralMatrix::getRow(int row)
 {
@@ -592,7 +602,7 @@ SVDDecomp::construct(const GeneralMatrix &A)
 }
 
 void
-SVDDecomp::solve(const GeneralMatrix &B, GeneralMatrix &X) const
+SVDDecomp::solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const
 {
   if (B.numRows() != U.numRows())
     throw SYLV_MES_EXCEPTION("Incompatible number of rows ");
diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh
index e0eb4c59f45385df552ad082de896c2506d36469..b88ac6413cbef6de0efe13d4dcae935ca29e9b78 100644
--- a/dynare++/sylv/cc/GeneralMatrix.hh
+++ b/dynare++/sylv/cc/GeneralMatrix.hh
@@ -6,6 +6,7 @@
 #define GENERAL_MATRIX_H
 
 #include "Vector.hh"
+#include "SylvException.hh"
 
 #include <algorithm>
 #include <memory>
@@ -17,7 +18,7 @@ class ConstGeneralMatrix
 {
   friend class GeneralMatrix;
 protected:
-  ConstVector data;
+  ConstVector data; // Has unit-stride
   int rows;
   int cols;
   int ld;
@@ -25,13 +26,22 @@ public:
   ConstGeneralMatrix(ConstVector d, int m, int n)
     : data(std::move(d)), rows(m), cols(n), ld(m)
   {
+    if (data.skip() > 1)
+      throw SYLV_MES_EXCEPTION("Vector must have unit-stride");
+    if (data.length() < m*n)
+      throw SYLV_MES_EXCEPTION("Vector is too small");
   }
+  ConstGeneralMatrix(const ConstGeneralMatrix &m) = default;
+  ConstGeneralMatrix(ConstGeneralMatrix &&m) = default;
   // Implicit conversion from ConstGeneralMatrix is ok, since it's cheap
   ConstGeneralMatrix(const GeneralMatrix &m);
   ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
   ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols);
   virtual ~ConstGeneralMatrix() = default;
 
+  ConstGeneralMatrix &operator=(const ConstGeneralMatrix &v) = delete;
+  ConstGeneralMatrix &operator=(ConstGeneralMatrix &&v) = delete;
+
   const double &
   get(int i, int j) const
   {
@@ -117,7 +127,7 @@ class GeneralMatrix
 {
   friend class ConstGeneralMatrix;
 protected:
-  Vector data;
+  Vector data; // Has unit-stride
   int rows;
   int cols;
   int ld;
@@ -126,17 +136,22 @@ public:
     : data(m*n), rows(m), cols(n), ld(m)
   {
   }
-  GeneralMatrix(const ConstVector &d, int m, int n)
-    : data(d), rows(m), cols(n), ld(m)
-  {
-  }
-  GeneralMatrix(Vector &d, int m, int n)
-    : data(d), rows(m), cols(n), ld(m)
+  GeneralMatrix(Vector d, int m, int n)
+    : data(std::move(d)), rows(m), cols(n), ld(m)
   {
+    if (data.skip() > 1)
+      throw SYLV_MES_EXCEPTION("Vector must have unit-stride");
+    if (data.length() < m*n)
+      throw SYLV_MES_EXCEPTION("Vector is too small");
   }
+
+  /* The copies will have ld==rows, for memory efficiency (hence we do not use
+     the default copy constructor) */
   GeneralMatrix(const GeneralMatrix &m);
   // We don't want implict conversion from ConstGeneralMatrix, since it's expensive
   explicit GeneralMatrix(const ConstGeneralMatrix &m);
+
+  GeneralMatrix(GeneralMatrix &&m) = default;
   GeneralMatrix(const GeneralMatrix &m, const char *dummy); // transpose
   GeneralMatrix(const ConstGeneralMatrix &m, const char *dummy); // transpose
   GeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
@@ -153,6 +168,8 @@ public:
 
   virtual ~GeneralMatrix() = default;
   GeneralMatrix &operator=(const GeneralMatrix &m) = default;
+  GeneralMatrix &operator=(GeneralMatrix &&m) = default;
+  GeneralMatrix &operator=(const ConstGeneralMatrix &m);
 
   const double &
   get(int i, int j) const
@@ -194,7 +211,7 @@ public:
   {
     return data;
   }
-  const Vector &
+  ConstVector
   getData() const
   {
     return data;
@@ -479,12 +496,12 @@ public:
   {
     return VT;
   }
-  void solve(const GeneralMatrix &B, GeneralMatrix &X) const;
+  void solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const;
   void
-  solve(const Vector &b, Vector &x) const
+  solve(const ConstVector &b, Vector &x) const
   {
     GeneralMatrix xmat(x, x.length(), 1);
-    solve(GeneralMatrix(b, b.length(), 1), xmat);
+    solve(ConstGeneralMatrix(b, b.length(), 1), xmat);
   }
 private:
   void construct(const GeneralMatrix &A);
diff --git a/dynare++/sylv/cc/GeneralSylvester.cc b/dynare++/sylv/cc/GeneralSylvester.cc
index 52ae68f1cd639e8c04d816b63700690ce9520c29..7cce950be42fc479fd5dd37b1cf346d793bf1aa4 100644
--- a/dynare++/sylv/cc/GeneralSylvester.cc
+++ b/dynare++/sylv/cc/GeneralSylvester.cc
@@ -15,8 +15,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
                                    const ConstVector &dc, const ConstVector &dd,
                                    const SylvParams &ps)
   : pars(ps),
-    mem_driver(pars, 1, m, n, ord), order(ord), a(da, n),
-    b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
+    mem_driver(pars, 1, m, n, ord), order(ord), a(Vector{da}, n),
+    b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(Vector{dd}, n, power(m, order)),
     solved(false)
 {
   init();
@@ -27,8 +27,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
                                    const ConstVector &dc, Vector &dd,
                                    const SylvParams &ps)
   : pars(ps),
-    mem_driver(pars, 0, m, n, ord), order(ord), a(da, n),
-    b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
+    mem_driver(pars, 0, m, n, ord), order(ord), a(Vector{da}, n),
+    b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(dd, n, power(m, order)),
     solved(false)
 {
   init();
@@ -39,8 +39,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
                                    const ConstVector &dc, const ConstVector &dd,
                                    bool alloc_for_check)
   : pars(alloc_for_check),
-    mem_driver(pars, 1, m, n, ord), order(ord), a(da, n),
-    b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
+    mem_driver(pars, 1, m, n, ord), order(ord), a(Vector{da}, n),
+    b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(Vector{dd}, n, power(m, order)),
     solved(false)
 {
   init();
@@ -51,8 +51,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
                                    const ConstVector &dc, Vector &dd,
                                    bool alloc_for_check)
   : pars(alloc_for_check),
-    mem_driver(pars, 0, m, n, ord), order(ord), a(da, n),
-    b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
+    mem_driver(pars, 0, m, n, ord), order(ord), a(Vector{da}, n),
+    b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(dd, n, power(m, order)),
     solved(false)
 {
   init();
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index f3a20d2fbc2bf22dbee2f825b47ee769c56d66fd..5f14e2d1b9397dd0fc5d76eb3b9fbe4f818def50 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -421,7 +421,7 @@ QuasiTriangular::QuasiTriangular(const QuasiTriangular &t)
 }
 
 QuasiTriangular::QuasiTriangular(const ConstVector &d, int d_size)
-  : SqSylvMatrix(d, d_size), diagonal(getData().base(), d_size)
+  : SqSylvMatrix(Vector{d}, d_size), diagonal(getData().base(), d_size)
 {
 }
 
diff --git a/dynare++/sylv/cc/QuasiTriangularZero.cc b/dynare++/sylv/cc/QuasiTriangularZero.cc
index fca2614b05dc6fef53abed23b517fdd422466a8c..260fc42ece1ed293688aa3fcecfd86b43679bf3e 100644
--- a/dynare++/sylv/cc/QuasiTriangularZero.cc
+++ b/dynare++/sylv/cc/QuasiTriangularZero.cc
@@ -8,15 +8,14 @@
 #include "SylvException.hh"
 
 #include <iostream>
-#include <utility>
 
 QuasiTriangularZero::QuasiTriangularZero(int num_zeros, const ConstVector &d,
                                          int d_size)
-  : QuasiTriangular(SqSylvMatrix(GeneralMatrix(d, num_zeros+d_size, d_size),
+  : QuasiTriangular(SqSylvMatrix(GeneralMatrix(Vector{d}, num_zeros+d_size, d_size),
                                  num_zeros, 0, d_size).getData(),
                     d_size),
     nz(num_zeros),
-    ru(GeneralMatrix(std::move(d), num_zeros+d_size, d_size), 0, 0, num_zeros, d_size)
+    ru(GeneralMatrix(Vector{d}, num_zeros+d_size, d_size), 0, 0, num_zeros, d_size)
 {
 }
 
diff --git a/dynare++/sylv/cc/SchurDecomp.hh b/dynare++/sylv/cc/SchurDecomp.hh
index f93e5aaf6e52aada6218611c274518450dfa8878..a8af976f5fbf4d403076692b337cd90d3ccbc8d4 100644
--- a/dynare++/sylv/cc/SchurDecomp.hh
+++ b/dynare++/sylv/cc/SchurDecomp.hh
@@ -50,7 +50,7 @@ class SchurDecompZero : public SchurDecomp
   GeneralMatrix ru; /* right upper matrix */
 public:
   SchurDecompZero(const GeneralMatrix &m);
-  const GeneralMatrix &
+  ConstGeneralMatrix
   getRU() const
   {
     return ru;
diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc
index 93e5807fec76810527ae9419584f85d26d473f17..ec16e76d7b38fcd88eec86321ad0b448d756cdc4 100644
--- a/dynare++/sylv/cc/SimilarityDecomp.cc
+++ b/dynare++/sylv/cc/SimilarityDecomp.cc
@@ -13,7 +13,7 @@
 
 SimilarityDecomp::SimilarityDecomp(const ConstVector &d, int d_size, double log10norm)
 {
-  SchurDecomp sd(SqSylvMatrix(d, d_size));
+  SchurDecomp sd(SqSylvMatrix(Vector{d}, d_size));
   q = std::make_unique<SqSylvMatrix>(sd.getQ());
   b = std::make_unique<BlockDiagonal>(sd.getT());
   invq = std::make_unique<SqSylvMatrix>(d_size);
diff --git a/dynare++/sylv/cc/SylvMatrix.hh b/dynare++/sylv/cc/SylvMatrix.hh
index b23ddd9e9de209610480dd71f80376ba0c1b994b..5e96a034a70d558a4d394225bb596fbd4c2c65e8 100644
--- a/dynare++/sylv/cc/SylvMatrix.hh
+++ b/dynare++/sylv/cc/SylvMatrix.hh
@@ -8,6 +8,8 @@
 #include "GeneralMatrix.hh"
 #include "KronVector.hh"
 
+#include <utility>
+
 class SqSylvMatrix;
 
 class SylvMatrix : public GeneralMatrix
@@ -17,12 +19,8 @@ public:
     : GeneralMatrix(m, n)
   {
   }
-  SylvMatrix(const ConstVector &d, int m, int n)
-    : GeneralMatrix(d, m, n)
-  {
-  }
-  SylvMatrix(Vector &d, int m, int n)
-    : GeneralMatrix(d, m, n)
+  SylvMatrix(Vector d, int m, int n)
+    : GeneralMatrix(std::move(d), m, n)
   {
   }
   SylvMatrix(const GeneralMatrix &m)
@@ -68,10 +66,7 @@ public:
   SqSylvMatrix(int m) : SylvMatrix(m, m)
   {
   }
-  SqSylvMatrix(const ConstVector &d, int m) : SylvMatrix(d, m, m)
-  {
-  }
-  SqSylvMatrix(Vector &d, int m) : SylvMatrix(d, m, m)
+  SqSylvMatrix(Vector d, int m) : SylvMatrix(std::move(d), m, m)
   {
   }
   SqSylvMatrix(const SqSylvMatrix &m) = default;
diff --git a/dynare++/sylv/cc/SymSchurDecomp.cc b/dynare++/sylv/cc/SymSchurDecomp.cc
index 0c7466e4c0df52ed50bf0acad760c921566ed990..c0f281a08583ac4330d0c52002d7fa75aba4aeef 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.cc
+++ b/dynare++/sylv/cc/SymSchurDecomp.cc
@@ -11,7 +11,7 @@
 #include <cmath>
 #include <vector>
 
-SymSchurDecomp::SymSchurDecomp(const GeneralMatrix &mata)
+SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata)
   : lambda(mata.numRows()), q(mata.numRows())
 {
   // check mata is square
diff --git a/dynare++/sylv/cc/SymSchurDecomp.hh b/dynare++/sylv/cc/SymSchurDecomp.hh
index 085649d24f84bd1f89f86590275b7d2452bbce62..b9d7e1ecccfc4d19fb02c1d96d22fd9ae882c4dc 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.hh
+++ b/dynare++/sylv/cc/SymSchurDecomp.hh
@@ -15,7 +15,7 @@ protected:
 public:
   /** Calculates A = Q*Lambda*Q^T, where A is assummed to be
    * symmetric and Lambda real diagonal, hence a vector. */
-  SymSchurDecomp(const GeneralMatrix &a);
+  SymSchurDecomp(const ConstGeneralMatrix &a);
   SymSchurDecomp(const SymSchurDecomp &ssd) = default;
   virtual ~SymSchurDecomp() = default;
   const Vector &
diff --git a/dynare++/sylv/testing/MMMatrix.cc b/dynare++/sylv/testing/MMMatrix.cc
index d9c42599030ef3eb51b45dda9fef1b389f9eba0d..9d3493d99cf47c0004e76e3266c5e554545b0ac5 100644
--- a/dynare++/sylv/testing/MMMatrix.cc
+++ b/dynare++/sylv/testing/MMMatrix.cc
@@ -24,7 +24,7 @@ MMMatrixIn::MMMatrixIn(const char *fname)
   if (2 != sscanf(buffer, "%d %d", &rows, &cols))
     throw MMException("Couldn't parse rows and cols\n");
   // read in data
-  data = std::shared_ptr<const double>(static_cast<double *>(operator new[](rows*cols*sizeof(double))), [](double *arr) { operator delete[](static_cast<void *>(arr)); });
+  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;
   int i = 0;
   while (fgets(buffer, 1000, fd) && i < len)
diff --git a/dynare++/sylv/testing/MMMatrix.hh b/dynare++/sylv/testing/MMMatrix.hh
index fc44ce0821dbbf7af82991be5662cbd544231dc0..3c6bb0e14ddc2ce79cbe500c89547a62f4a6807c 100644
--- a/dynare++/sylv/testing/MMMatrix.hh
+++ b/dynare++/sylv/testing/MMMatrix.hh
@@ -33,16 +33,16 @@ public:
 
 class MMMatrixIn : public MallocAllocator
 {
-  std::shared_ptr<const double> data;
+  std::shared_ptr<double> data;
   int rows;
   int cols;
 public:
   MMMatrixIn(const char *fname);
   ~MMMatrixIn() = default;
-  ConstVector
+  Vector
   getData() const
   {
-    return ConstVector{data, size()};
+    return Vector{data, size()};
   }
   int
   size() const
diff --git a/dynare++/tl/cc/kron_prod.cc b/dynare++/tl/cc/kron_prod.cc
index d15f746a56cd4f28c26919a49c067e30ec29781d..6191cccdbe6ee6c4f3f6663e96fe21e8cb1704db 100644
--- a/dynare++/tl/cc/kron_prod.cc
+++ b/dynare++/tl/cc/kron_prod.cc
@@ -234,7 +234,7 @@ KronProdIAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
 }
 
 /* Here we multiply $B\cdot(A_1\otimes\ldots\otimes A_n)$. First we
-   multiply $B\cdot(A_1\otimes)$, then this is multiplied by all
+   multiply $B\cdot(A_1\otimes I)$, then this is multiplied by all
    $I\otimes A_i\otimes I$, and finally by $I\otimes A_n$.
 
    If the dimension of the Kronecker product is only 1, then we multiply
@@ -293,7 +293,7 @@ KronProdAll::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
     }
   else
     {
-      last = new TwoDMatrix(in.nrows(), in.ncols(), in.getData());
+      last = new TwoDMatrix(in.nrows(), in.ncols(), Vector{in.getData()});
     }
 
   // perform intermediate multiplications IAI
diff --git a/dynare++/tl/cc/twod_matrix.cc b/dynare++/tl/cc/twod_matrix.cc
index a88896ff5cbc34bc59d0548487fed812758f7c0f..c71393f4415f7411b7848a4197ca681bee710b56 100644
--- a/dynare++/tl/cc/twod_matrix.cc
+++ b/dynare++/tl/cc/twod_matrix.cc
@@ -54,6 +54,13 @@ ConstTwoDMatrix::writeMat(mat_t *fd, const char *vname) const
   delete[] data;
 }
 
+TwoDMatrix &
+TwoDMatrix::operator=(const ConstTwoDMatrix &m)
+{
+  GeneralMatrix::operator=(m);
+  return *this;
+}
+
 void
 TwoDMatrix::copyRow(int from, int to)
 {
diff --git a/dynare++/tl/cc/twod_matrix.hh b/dynare++/tl/cc/twod_matrix.hh
index 43c0f3bc893a93f960bd8796a463b8154a8a5889..d0d463f958b20690bd7b8c7554963a14ad179fe7 100644
--- a/dynare++/tl/cc/twod_matrix.hh
+++ b/dynare++/tl/cc/twod_matrix.hh
@@ -34,6 +34,8 @@ public:
     : ConstGeneralMatrix(std::move(d), m, n)
   {
   }
+  ConstTwoDMatrix(const ConstTwoDMatrix &m) = default;
+  ConstTwoDMatrix(ConstTwoDMatrix &&m) = default;
   // Implicit conversion from TwoDMatrix is ok, since it's cheap
   ConstTwoDMatrix(const TwoDMatrix &m);
   ConstTwoDMatrix(const TwoDMatrix &m, int first_col, int num);
@@ -47,6 +49,9 @@ public:
   ~ConstTwoDMatrix()
   override = default;
 
+  ConstTwoDMatrix &operator=(const ConstTwoDMatrix &v) = delete;
+  ConstTwoDMatrix &operator=(ConstTwoDMatrix &&v) = delete;
+
   int
   nrows() const
   {
@@ -69,16 +74,14 @@ public:
 class TwoDMatrix : public GeneralMatrix
 {
 public:
+  TwoDMatrix(const TwoDMatrix &m) = default;
+  TwoDMatrix(TwoDMatrix &&m) = default;
   TwoDMatrix(int r, int c)
     : GeneralMatrix(r, c)
   {
   }
-  TwoDMatrix(int r, int c, Vector &d)
-    : GeneralMatrix(d, r, c)
-  {
-  }
-  TwoDMatrix(int r, int c, const ConstVector &d)
-    : GeneralMatrix(d, r, c)
+  TwoDMatrix(int r, int c, Vector d)
+    : GeneralMatrix(std::move(d), r, c)
   {
   }
   TwoDMatrix(const GeneralMatrix &m)
@@ -122,8 +125,11 @@ public:
     : GeneralMatrix(a, b)
   {
   }
-  ~TwoDMatrix()
-  override = default;
+  ~TwoDMatrix() override = default;
+
+  TwoDMatrix &operator=(const TwoDMatrix &m) = default;
+  TwoDMatrix &operator=(TwoDMatrix &&m) = default;
+  TwoDMatrix &operator=(const ConstTwoDMatrix &m);
 
   int
   nrows() const
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index c59fe8325c0d1ca49b9a1121a1055567d4e415b5..4d5e5c25f338fa828b8068e2ed9fde85531a46e0 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -200,19 +200,19 @@ extern "C" {
         const mxArray *g1 = prhs[3];
         int m = (int) mxGetM(g1);
         int n = (int) mxGetN(g1);
-        g1m = new TwoDMatrix(m, n, ConstVector{g1});
+        g1m = new TwoDMatrix(m, n, Vector{ConstVector{g1}});
         if (nrhs > 4)
           {
             const mxArray *g2 = prhs[4];
             int m = (int) mxGetM(g2);
             int n = (int) mxGetN(g2);
-            g2m = new TwoDMatrix(m, n, ConstVector{g2});
+            g2m = new TwoDMatrix(m, n, Vector{ConstVector{g2}});
             if (nrhs > 5)
               {
                 const mxArray *g3 = prhs[5];
                 int m = (int) mxGetM(g3);
                 int n = (int) mxGetN(g3);
-                g3m = new TwoDMatrix(m, n, ConstVector{g3});
+                g3m = new TwoDMatrix(m, n, Vector{ConstVector{g3}});
               }
           }
       }