diff --git a/dynare++/integ/cc/vector_function.cc b/dynare++/integ/cc/vector_function.cc
index 92afd6ee7653f58afa303664cea77413fee4c9d0..d93881a15221a2962d94e08e4f951161a63ce6ee 100644
--- a/dynare++/integ/cc/vector_function.cc
+++ b/dynare++/integ/cc/vector_function.cc
@@ -120,12 +120,12 @@ GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
 {
   A = vcov;
 
-  lapack_int rows = A.numRows();
+  lapack_int rows = A.numRows(), lda = A.getLD();
   for (int i = 0; i < rows; i++)
     for (int j = i+1; j < rows; j++)
       A.get(i, j) = 0.0;
 
   lapack_int info;
-  dpotrf("L", &rows, A.base(), &rows, &info);
+  dpotrf("L", &rows, A.base(), &lda, &info);
   // todo: raise if |info!=1|
 }
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index e1998d0d77bb5165ae64ab6a61a0100762e5092f..f61d937b5a21ca2e7eafc2947f055f3ac95ad337 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -563,13 +563,13 @@ void
 RandomShockRealization::choleskyFactor(const TwoDMatrix &v)
 {
   factor = v;
-  lapack_int rows = factor.nrows();
+  lapack_int rows = factor.nrows(), lda = factor.getLD();
   for (int i = 0; i < rows; i++)
     for (int j = i+1; j < rows; j++)
       factor.get(i, j) = 0.0;
   lapack_int info;
 
-  dpotrf("L", &rows, factor.base(), &rows, &info);
+  dpotrf("L", &rows, factor.base(), &lda, &info);
   KORD_RAISE_IF(info != 0,
                 "Info!=0 in RandomShockRealization::choleskyFactor");
 }
diff --git a/dynare++/kord/first_order.cc b/dynare++/kord/first_order.cc
index 7986481337f058d91f1fb3de205d56788cce3881..733d51670d31d98c7cc35fc2eddeb0e8d40682c4 100644
--- a/dynare++/kord/first_order.cc
+++ b/dynare++/kord/first_order.cc
@@ -145,6 +145,7 @@ FirstOrder::solve(const TwoDMatrix &fd)
   matD.place(fyplus, 0, ypart.nys()+ypart.nstat);
   for (int i = 0; i < ypart.nboth; i++)
     matD.get(ypart.ny()+i, ypart.npred+i) = 1.0;
+  lapack_int ldb = matD.getLD();
 
   // form matrix $E$
   TwoDMatrix matE(n, n);
@@ -155,18 +156,20 @@ FirstOrder::solve(const TwoDMatrix &fd)
   for (int i = 0; i < ypart.nboth; i++)
     matE.get(ypart.ny()+i, ypart.nys()+ypart.nstat+i) = -1.0;
   matE.mult(-1.0);
+  lapack_int lda = matE.getLD();
 
   // solve generalized Schur
   TwoDMatrix vsl(n, n);
   TwoDMatrix vsr(n, n);
+  lapack_int ldvsl = vsl.getLD(), ldvsr = vsr.getLD();
   lapack_int lwork = 100*n+16;
   Vector work(lwork);
   auto *bwork = new lapack_int[n];
   lapack_int info;
   lapack_int sdim2 = sdim;
-  dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &n,
-        matD.getData().base(), &n, &sdim2, alphar.base(), alphai.base(),
-        beta.base(), vsl.getData().base(), &n, vsr.getData().base(), &n,
+  dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &lda,
+        matD.getData().base(), &ldb, &sdim2, alphar.base(), alphai.base(),
+        beta.base(), vsl.getData().base(), &ldvsl, vsr.getData().base(), &ldvsr,
         work.base(), &lwork, bwork, &info);
   if (info)
     {
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index c9bcef5be6fa26944dd5341eb16388c6a5ff7659..89893f7bcfb46c3bd72e570e01e0e6a072998c65 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -17,9 +17,9 @@ void
 PLUMatrix::calcPLU()
 {
   lapack_int info;
-  lapack_int rows = nrows();
+  lapack_int rows = nrows(), lda = ld;
   inv = (const Vector &) getData();
-  dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
+  dgetrf(&rows, &rows, inv.base(), &lda, ipiv, &info);
 }
 
 /* Here we just call the LAPACK machinery to multiply by the inverse. */
@@ -30,11 +30,13 @@ PLUMatrix::multInv(TwoDMatrix &m) const
   KORD_RAISE_IF(m.nrows() != ncols(),
                 "The matrix is not square in PLUMatrix::multInv");
   lapack_int info;
+  lapack_int lda = ld;
   lapack_int mcols = m.ncols();
   lapack_int mrows = m.nrows();
+  lapack_int ldb = m.getLD();
   double *mbase = m.getData().base();
-  dgetrs("N", &mrows, &mcols, inv.base(), &mrows, ipiv,
-         mbase, &mrows, &info);
+  dgetrs("N", &mrows, &mcols, inv.base(), &lda, ipiv,
+         mbase, &ldb, &info);
   KORD_RAISE_IF(info != 0,
                 "Info!=0 in PLUMatrix::multInv");
 }
diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc
index 7a659c97348ec329b45f8f15a22d2699c8adfdd0..7ba83d783ae6a0c0749918238a1a4da5cb66cbc3 100644
--- a/dynare++/sylv/cc/GeneralMatrix.cc
+++ b/dynare++/sylv/cc/GeneralMatrix.cc
@@ -448,7 +448,7 @@ ConstGeneralMatrix::multVecTrans(double a, Vector &x, double b,
       blas_int mm = rows;
       blas_int nn = cols;
       double alpha = b;
-      blas_int lda = rows;
+      blas_int lda = ld;
       blas_int incx = d.skip();
       double beta = a;
       blas_int incy = x.skip();
@@ -471,9 +471,9 @@ ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld
       GeneralMatrix inv(*this);
       std::vector<lapack_int> ipiv(rows);
       lapack_int info;
-      lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld;
-      dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv.data(), &info);
-      dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv.data(), d,
+      lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld, lda = inv.ld;
+      dgetrf(&rows2, &rows2, inv.getData().base(), &lda, ipiv.data(), &info);
+      dgetrs(trans, &rows2, &mcols2, inv.base(), &lda, ipiv.data(), d,
              &mld2, &info);
     }
 }
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index b0fda9f97090f18b7386d2fc94efa1e6e83b7f29..f3a20d2fbc2bf22dbee2f825b47ee769c56d66fd 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -432,8 +432,9 @@ QuasiTriangular::QuasiTriangular(int p, const QuasiTriangular &t)
   blas_int d_size = diagonal.getSize();
   double alpha = 1.0;
   double beta = 0.0;
+  blas_int lda = t.getLD(), ldb = t.getLD(), ldc = ld;
   dgemm("N", "N", &d_size, &d_size, &d_size, &alpha, aux.base(),
-        &d_size, t.getData().base(), &d_size, &beta, getData().base(), &d_size);
+        &lda, t.getData().base(), &ldb, &beta, getData().base(), &ldc);
 }
 
 QuasiTriangular::QuasiTriangular(const SchurDecomp &decomp)
@@ -590,7 +591,7 @@ QuasiTriangular::solvePre(Vector &x, double &eig_min)
     }
 
   blas_int nn = diagonal.getSize();
-  blas_int lda = diagonal.getSize();
+  blas_int lda = ld;
   blas_int incx = x.skip();
   dtrsv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
 }
@@ -616,7 +617,7 @@ QuasiTriangular::solvePreTrans(Vector &x, double &eig_min)
     }
 
   blas_int nn = diagonal.getSize();
-  blas_int lda = diagonal.getSize();
+  blas_int lda = ld;
   blas_int incx = x.skip();
   dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
 }
@@ -627,7 +628,7 @@ QuasiTriangular::multVec(Vector &x, const ConstVector &b) const
 {
   x = b;
   blas_int nn = diagonal.getSize();
-  blas_int lda = diagonal.getSize();
+  blas_int lda = ld;
   blas_int incx = x.skip();
   dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
   for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
@@ -645,7 +646,7 @@ QuasiTriangular::multVecTrans(Vector &x, const ConstVector &b) const
 {
   x = b;
   blas_int nn = diagonal.getSize();
-  blas_int lda = diagonal.getSize();
+  blas_int lda = ld;
   blas_int incx = x.skip();
   dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
   for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
diff --git a/dynare++/sylv/cc/SchurDecomp.cc b/dynare++/sylv/cc/SchurDecomp.cc
index 7a101d47d43f1e9046f377ed017e43d9b3164f6a..8cc1abc751effca21a2dcbda81be53b50246a312 100644
--- a/dynare++/sylv/cc/SchurDecomp.cc
+++ b/dynare++/sylv/cc/SchurDecomp.cc
@@ -13,13 +13,14 @@ SchurDecomp::SchurDecomp(const SqSylvMatrix &m)
 {
   lapack_int rows = m.numRows();
   SqSylvMatrix auxt(m);
+  lapack_int lda = auxt.getLD(), ldvs = q.getLD();
   lapack_int sdim;
   std::vector<double> wr(rows), wi(rows);
   lapack_int lwork = 6*rows;
   std::vector<double> work(lwork);
   lapack_int info;
-  dgees("V", "N", nullptr, &rows, auxt.base(), &rows, &sdim,
-        wr.data(), wi.data(), q.base(), &rows,
+  dgees("V", "N", nullptr, &rows, auxt.base(), &lda, &sdim,
+        wr.data(), wi.data(), q.base(), &ldvs,
         work.data(), &lwork, nullptr, &info);
   t_storage = std::make_unique<QuasiTriangular>(auxt.getData(), rows);
   t = t_storage.get();
diff --git a/dynare++/sylv/cc/SchurDecompEig.cc b/dynare++/sylv/cc/SchurDecompEig.cc
index 3c5955d32625d01005abbfe41cd971414c6c7b6b..5193ed7bf127df2ed495923dde50c727d2d0718f 100644
--- a/dynare++/sylv/cc/SchurDecompEig.cc
+++ b/dynare++/sylv/cc/SchurDecompEig.cc
@@ -54,12 +54,12 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
   itadd = it;
   --itadd;
 
-  lapack_int n = getDim();
+  lapack_int n = getDim(), ldt = getT().getLD(), ldq = getQ().getLD();
   lapack_int ifst = (*it).getIndex() + 1;
   lapack_int ilst = (*itadd).getIndex() + 1;
   std::vector<double> work(n);
   lapack_int info;
-  dtrexc("V", &n, getT().base(), &n, getQ().base(), &n, &ifst, &ilst, work.data(),
+  dtrexc("V", &n, getT().base(), &ldt, getQ().base(), &ldq, &ifst, &ilst, work.data(),
          &info);
   if (info < 0)
     throw SYLV_MES_EXCEPTION("Wrong argument to dtrexc.");
diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc
index 0da2ebf7e8c4584949451092d1972557cf3f9471..93e5807fec76810527ae9419584f85d26d473f17 100644
--- a/dynare++/sylv/cc/SimilarityDecomp.cc
+++ b/dynare++/sylv/cc/SimilarityDecomp.cc
@@ -51,9 +51,10 @@ SimilarityDecomp::solveX(diag_iter start, diag_iter end,
   lapack_int isgn = -1;
   lapack_int m = A.numRows();
   lapack_int n = B.numRows();
+  lapack_int lda = A.getLD(), ldb = B.getLD();
   double scale;
   lapack_int info;
-  dtrsyl("N", "N", &isgn, &m, &n, A.base(), &m, B.base(), &n,
+  dtrsyl("N", "N", &isgn, &m, &n, A.base(), &lda, B.base(), &ldb,
          C.base(), &m, &scale, &info);
   if (info < -1)
     throw SYLV_MES_EXCEPTION("Wrong parameter to LAPACK dtrsyl.");
diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc
index 6af0480c03dacfc3663d8208580aa8d65d664944..1e0b3241bb1c80b16b3e53f30e2a6a3acd5db0f8 100644
--- a/dynare++/sylv/cc/SylvMatrix.cc
+++ b/dynare++/sylv/cc/SylvMatrix.cc
@@ -233,27 +233,27 @@ SqSylvMatrix::multInvLeft2(GeneralMatrix &a, GeneralMatrix &b,
   Vector inv(data);
   std::vector<lapack_int> ipiv(rows);
   lapack_int info;
-  lapack_int rows2 = rows;
-  dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv.data(), &info);
+  lapack_int rows2 = rows, lda = ld;
+  dgetrf(&rows2, &rows2, inv.base(), &lda, ipiv.data(), &info);
   // solve a
   lapack_int acols = a.numCols();
   double *abase = a.base();
-  dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv.data(),
+  dgetrs("N", &rows2, &acols, inv.base(), &lda, ipiv.data(),
          abase, &rows2, &info);
   // solve b
   lapack_int bcols = b.numCols();
   double *bbase = b.base();
-  dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv.data(),
+  dgetrs("N", &rows2, &bcols, inv.base(), &lda, ipiv.data(),
          bbase, &rows2, &info);
 
   // condition numbers
   std::vector<double> work(4*rows);
   std::vector<lapack_int> iwork(rows);
   double norm1 = getNorm1();
-  dgecon("1", &rows2, inv.base(), &rows2, &norm1, &rcond1,
+  dgecon("1", &rows2, inv.base(), &lda, &norm1, &rcond1,
          work.data(), iwork.data(), &info);
   double norminf = getNormInf();
-  dgecon("I", &rows2, inv.base(), &rows2, &norminf, &rcondinf,
+  dgecon("I", &rows2, inv.base(), &lda, &norminf, &rcondinf,
          work.data(), iwork.data(), &info);
 }