k-order DLL: various simplifications

mex/build/k_order_perturbation.am

@@ -15,7 +15,6 @@ nodist_k_order_perturbation_SOURCES = \
 	$(TOPDIR)/k_ord_dynare.hh \
 	$(TOPDIR)/dynamic_dll.cc \
 	$(TOPDIR)/dynamic_dll.hh \
-	$(TOPDIR)/dynamic_abstract_class.cc \
 	$(TOPDIR)/dynamic_abstract_class.hh \
 	$(TOPDIR)/dynamic_m.cc \
 	$(TOPDIR)/dynamic_m.hh

mex/sources/k_order_perturbation/dynamic_abstract_class.cc

@@ -1,72 +0,0 @@
-/*
- * Copyright (C) 2010-2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include <assert.h>
-
-#include "dynamic_abstract_class.hh"
-
-void
-DynamicModelAC::copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols)
-{
-  assert(rows == tdm->nrows());
-  assert(cols == tdm->ncols());
-
-  int dmIdx = 0;
-  for (int j = 0; j < cols; j++)
-    for (int i = 0; i < rows; i++)
-      tdm->get(i, j) = dm[dmIdx++];
-}
-
-void
-DynamicModelAC::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm)
-{
-  int totalCols = mxGetN(sparseMat);
-  mwIndex *rowIdxVector = mxGetIr(sparseMat);
-  mwSize sizeRowIdxVector = mxGetNzmax(sparseMat);
-  mwIndex *colIdxVector = mxGetJc(sparseMat);
-
-  assert(tdm->ncols() == 3);
-  assert(tdm->nrows() == sizeRowIdxVector);
-
-  double *ptr = mxGetPr(sparseMat);
-
-  int rind = 0;
-  int output_row = 0;
-
-  for (int i = 0; i < totalCols; i++)
-    for (int j = 0; j < static_cast<int>((colIdxVector[i+1]-colIdxVector[i])); j++, rind++)
-      {
-        tdm->get(output_row, 0) = rowIdxVector[rind] + 1;
-        tdm->get(output_row, 1) = i + 1;
-        tdm->get(output_row, 2) = ptr[rind];
-        output_row++;
-      }
-
-  /* If there are less elements than Nzmax (that might happen if some
-     derivative is symbolically not zero but numerically zero at the evaluation
-     point), then fill in the matrix with empty entries, that will be
-     recognized as such by KordpDynare::populateDerivativesContainer() */
-  while (output_row < static_cast<int>(sizeRowIdxVector))
-    {
-      tdm->get(output_row, 0) = 0;
-      tdm->get(output_row, 1) = 0;
-      tdm->get(output_row, 2) = 0;
-      output_row++;
-    }
-}

mex/sources/k_order_perturbation/dynamic_abstract_class.hh

@@ -20,15 +20,13 @@
 #ifndef _DYNAMICMODELAC_HH
 #define _DYNAMICMODELAC_HH
 
-#include "k_ord_dynare.hh"
+#include "twod_matrix.hh"
 
 class DynamicModelAC
 {
 public:
   virtual ~DynamicModelAC() = default;
-  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm);
-  static void copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols);
   virtual void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-                    Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false) = 0;
+                    Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) = 0;
 };
 #endif

mex/sources/k_order_perturbation/dynamic_dll.cc

@@ -19,6 +19,10 @@
 
 #include "dynamic_dll.hh"
 
+#include "dynare_exception.hh"
+
+#include <iostream>
+
 DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
 {
   std::string fName;
@@ -27,71 +31,51 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName) noexcept(false)
 #endif
   fName += "+" + modName + "/dynamic" + MEXEXT;
 
-  try
-    {
 #if defined(__CYGWIN32__) || defined(_WIN32)
-      dynamicHinstance = LoadLibrary(fName.c_str());
+  dynamicHinstance = LoadLibrary(fName.c_str());
-      if (!dynamicHinstance)
+  if (!dynamicHinstance)
-        throw 1;
+    throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": can't dynamically load the file");
-      ntt = reinterpret_cast<int *>(GetProcAddress(dynamicHinstance, "ntt"));
+  ntt = reinterpret_cast<int *>(GetProcAddress(dynamicHinstance, "ntt"));
-      dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid_tt"));
+  dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid_tt"));
-      dynamic_resid = reinterpret_cast<dynamic_resid_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid"));
+  dynamic_resid = reinterpret_cast<dynamic_resid_fct>(GetProcAddress(dynamicHinstance, "dynamic_resid"));
-      dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1_tt"));
+  dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1_tt"));
-      dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1"));
+  dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(GetProcAddress(dynamicHinstance, "dynamic_g1"));
-      dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2_tt"));
+  dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2_tt"));
-      dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2"));
+  dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(GetProcAddress(dynamicHinstance, "dynamic_g2"));
-      dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3_tt"));
+  dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3_tt"));
-      dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3"));
+  dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(GetProcAddress(dynamicHinstance, "dynamic_g3"));
-      if (!ntt || !dynamic_resid_tt || !dynamic_resid
+  if (!ntt || !dynamic_resid_tt || !dynamic_resid
-          || !dynamic_g1_tt || !dynamic_g1
+      || !dynamic_g1_tt || !dynamic_g1
-          || !dynamic_g2_tt || !dynamic_g2
+      || !dynamic_g2_tt || !dynamic_g2
-          || !dynamic_g3_tt || !dynamic_g3)
+      || !dynamic_g3_tt || !dynamic_g3)
-        {
+    {
-          FreeLibrary(dynamicHinstance); // Free the library
+      FreeLibrary(dynamicHinstance); // Free the library
-          throw 2;
+      throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": can't locate the relevant dynamic symbols within the MEX file");
-        }
+    }
 #else // Linux or Mac
-      dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW);
+  dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW);
-      if (dynamicHinstance == nullptr)
+  if (!dynamicHinstance)
-        {
+    throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": " + dlerror());
-          std::cerr << dlerror() << std::endl;
+  ntt = reinterpret_cast<int *>(dlsym(dynamicHinstance, "ntt"));
-          throw 1;
+  dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_resid_tt"));
-        }
+  dynamic_resid = reinterpret_cast<dynamic_resid_fct>(dlsym(dynamicHinstance, "dynamic_resid"));
-      ntt = reinterpret_cast<int *>(dlsym(dynamicHinstance, "ntt"));
+  dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g1_tt"));
-      dynamic_resid_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_resid_tt"));
+  dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(dlsym(dynamicHinstance, "dynamic_g1"));
-      dynamic_resid = reinterpret_cast<dynamic_resid_fct>(dlsym(dynamicHinstance, "dynamic_resid"));
+  dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g2_tt"));
-      dynamic_g1_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g1_tt"));
+  dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(dlsym(dynamicHinstance, "dynamic_g2"));
-      dynamic_g1 = reinterpret_cast<dynamic_g1_fct>(dlsym(dynamicHinstance, "dynamic_g1"));
+  dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g3_tt"));
-      dynamic_g2_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g2_tt"));
+  dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(dlsym(dynamicHinstance, "dynamic_g3"));
-      dynamic_g2 = reinterpret_cast<dynamic_g2_fct>(dlsym(dynamicHinstance, "dynamic_g2"));
+  if (!ntt || !dynamic_resid_tt || !dynamic_resid
-      dynamic_g3_tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, "dynamic_g3_tt"));
+      || !dynamic_g1_tt || !dynamic_g1
-      dynamic_g3 = reinterpret_cast<dynamic_g3_fct>(dlsym(dynamicHinstance, "dynamic_g3"));
+      || !dynamic_g2_tt || !dynamic_g2
-      if (!ntt || !dynamic_resid_tt || !dynamic_resid
+      || !dynamic_g3_tt || !dynamic_g3)
-          || !dynamic_g1_tt || !dynamic_g1
+    {
-          || !dynamic_g2_tt || !dynamic_g2
+      dlclose(dynamicHinstance); // Free the library
-          || !dynamic_g3_tt || !dynamic_g3)
+      throw DynareException(__FILE__, __LINE__, "Error when loading " + fName + ": " + dlerror());
-        {
+    }
-          dlclose(dynamicHinstance); // Free the library
-          std::cerr << dlerror() << std::endl;
-          throw 2;
-        }
 #endif
 
-    }
+  tt = std::make_unique<double[]>(*ntt);
-  catch (int i)
-    {
-      std::string msg{"Error when loading " + fName + " ("};
-      if (i == 1)
-        msg += "can't dynamically load the file";
-      if (i == 2)
-        msg += "can't locate the relevant dynamic symbols within the MEX file";
-      msg += ")";
-      throw DynareException(__FILE__, __LINE__, msg);
-    }
-  catch (...)
-    {
-      throw DynareException(__FILE__, __LINE__, "Can't find the relevant dynamic symbols in " + fName);
-    }
 }
 
 DynamicModelDLL::~DynamicModelDLL()
@@ -112,21 +96,19 @@ void
 DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
                       Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false)
 {
-  double *T = new double[*ntt];
+  dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-  dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
+  dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), residual.base());
-  dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, residual.base());
   if (g1 || g2 || g3)
-    dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
+    dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
   if (g1)
-    dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g1->base());
+    dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g1->base());
   if (g2 || g3)
-    dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
+    dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
   if (g2)
-    dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g2->base());
+    dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g2->base());
   if (g3)
     {
-      dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T);
+      dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-      dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, T, g3->base());
+      dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g3->base());
     }
-  delete[] T;
 }

mex/sources/k_order_perturbation/dynamic_dll.hh

@@ -30,9 +30,9 @@
 #endif
 
 #include <string>
+#include <memory>
 
 #include "dynamic_abstract_class.hh"
-#include "dynare_exception.hh"
 
 using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);
 using dynamic_resid_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *residual);
@@ -58,13 +58,14 @@ private:
 #else
   void *dynamicHinstance; // and in Linux or Mac
 #endif
+  std::unique_ptr<double[]> tt; // Vector of temporary terms
 
 public:
   // construct and load Dynamic model DLL
-  explicit DynamicModelDLL(const std::string &fname) noexcept(false);
+  explicit DynamicModelDLL(const std::string &fname);
   virtual ~DynamicModelDLL();
 
   void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false);
+            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
 };
 #endif

mex/sources/k_order_perturbation/dynamic_m.cc

@@ -18,6 +18,9 @@
  */
 
 #include <algorithm>
+#include <cassert>
+
+#include "dynare_exception.hh"
 
 #include "dynamic_m.hh"
 
@@ -26,10 +29,49 @@ DynamicModelMFile::DynamicModelMFile(const std::string &modName) noexcept(false)
 {
 }
 
+void
+DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm)
+{
+  int totalCols = mxGetN(sparseMat);
+  mwIndex *rowIdxVector = mxGetIr(sparseMat);
+  mwSize sizeRowIdxVector = mxGetNzmax(sparseMat);
+  mwIndex *colIdxVector = mxGetJc(sparseMat);
+
+  assert(tdm->ncols() == 3);
+  assert(tdm->nrows() == sizeRowIdxVector);
+
+  double *ptr = mxGetPr(sparseMat);
+
+  int rind = 0;
+  int output_row = 0;
+
+  for (int i = 0; i < totalCols; i++)
+    for (int j = 0; j < static_cast<int>((colIdxVector[i+1]-colIdxVector[i])); j++, rind++)
+      {
+        tdm->get(output_row, 0) = rowIdxVector[rind] + 1;
+        tdm->get(output_row, 1) = i + 1;
+        tdm->get(output_row, 2) = ptr[rind];
+        output_row++;
+      }
+
+  /* If there are less elements than Nzmax (that might happen if some
+     derivative is symbolically not zero but numerically zero at the evaluation
+     point), then fill in the matrix with empty entries, that will be
+     recognized as such by KordpDynare::populateDerivativesContainer() */
+  while (output_row < static_cast<int>(sizeRowIdxVector))
+    {
+      tdm->get(output_row, 0) = 0;
+      tdm->get(output_row, 1) = 0;
+      tdm->get(output_row, 2) = 0;
+      output_row++;
+    }
+}
+
 void
 DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
                         Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false)
 {
+  constexpr int nlhs_dynamic = 4, nrhs_dynamic = 5;
   mxArray *prhs[nrhs_dynamic], *plhs[nlhs_dynamic];
 
   prhs[0] = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
@@ -48,8 +90,11 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
     throw DynareException(__FILE__, __LINE__, "Trouble calling " + DynamicMFilename);
 
   residual = Vector{plhs[0]};
-  copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, static_cast<int>(mxGetM(plhs[1])),
+
-                            static_cast<int>(mxGetN(plhs[1])));
+  assert(static_cast<int>(mxGetM(plhs[1])) == g1->nrows());
+  assert(static_cast<int>(mxGetN(plhs[1])) == g1->ncols());
+  std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), g1->base());
+
   if (g2)
     unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
   if (g3)

mex/sources/k_order_perturbation/dynamic_m.hh

@@ -33,12 +33,11 @@ class DynamicModelMFile : public DynamicModelAC
 {
 private:
   const std::string DynamicMFilename;
-  constexpr static int nlhs_dynamic = 4;
+  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm);
-  constexpr static int nrhs_dynamic = 5;
 public:
-  explicit DynamicModelMFile(const std::string &modName) noexcept(false);
+  explicit DynamicModelMFile(const std::string &modName);
   virtual ~DynamicModelMFile() = default;
   void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false);
+            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
 };
 #endif

mex/sources/k_order_perturbation/k_ord_dynare.cc

@@ -19,74 +19,43 @@
 
 // GP, based on work by O.Kamenik
 
-#include <vector>
+#include "k_ord_dynare.hh"
-#include "first_order.hh"
 #include "dynamic_abstract_class.hh"
+#include "dynare_exception.hh"
 
-#include <cmath>
+#include <utility>
-#include <sstream>
 
-#include <iostream>
-#include <fstream>
 
-/**************************************************************************************/
+KordpDynare::KordpDynare(const std::vector<std::string> &endo,
-/*       Dynare DynamicModel class                                                                 */
-/**************************************************************************************/
-
-KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
                          const std::vector<std::string> &exo, int nexog, int npar,
                          Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
-                         int npred, int nforw, int nboth, const int jcols, const Vector &nnzd,
+                         int npred, int nforw, int nboth, const Vector &nnzd,
-                         const int nsteps, int norder,
+                         int nsteps, int norder,
-                         Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
+                         Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
-                         const std::vector<int> &var_order, const TwoDMatrix &llincidence, double criterium) noexcept(false) :
+                         const std::vector<int> &var_order, const TwoDMatrix &llincidence,
-  nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
-  nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
-  nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
-  md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
-  ll_Incidence{llincidence}, qz_criterium{criterium},
-  dynamicModelFile{move(dynamicModelFile_arg)}
-{
-  ReorderDynareJacobianIndices();
-
-  //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
-  for (int iord = 1; iord <= nOrder; iord++)
-    md.insert(std::make_unique<FSSparseTensor>(iord, nY+nYs+nYss+nExog, nY));
-}
-
-KordpDynare::KordpDynare(const std::vector<std::string> &endo, int num_endo,
-                         const std::vector<std::string> &exo, int nexog, int npar,
-                         Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
-                         int npred, int nforw, int nboth, const int jcols, const Vector &nnzd,
-                         const int nsteps, int norder,
-                         Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
-                         const std::vector<int> &var_order, const TwoDMatrix &llincidence, double criterium,
                          std::unique_ptr<TwoDMatrix> g1_arg, std::unique_ptr<TwoDMatrix> g2_arg,
-                         std::unique_ptr<TwoDMatrix> g3_arg) noexcept(false) :
+                         std::unique_ptr<TwoDMatrix> g3_arg) :
   nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
-  nYs{npred + nboth}, nYss{nboth + nforw}, nY{num_endo}, nJcols{jcols}, NNZD{nnzd}, nSteps{nsteps},
+  nYs{npred + nboth}, nYss{nboth + nforw}, nY{nstat + npred + nboth + nforw},
+  nJcols{nExog+nY+nYs+nYss}, NNZD{nnzd}, nSteps{nsteps},
   nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
-  md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl}, ss_tol{sstol}, varOrder{var_order},
+  md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl},
-  ll_Incidence{llincidence}, qz_criterium{criterium},
+  ll_Incidence{llincidence},
   g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)},
-  dynamicModelFile{move(dynamicModelFile_arg)}
+  dynamicModelFile{std::move(dynamicModelFile_arg)}
 {
-  ReorderDynareJacobianIndices();
+  computeJacobianPermutation(var_order);
-
-  //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
-  for (int iord = 1; iord <= nOrder; iord++)
-    md.insert(std::make_unique<FSSparseTensor>(iord, nY+nYs+nYss+nExog, nY));
 }
 
 void
 KordpDynare::solveDeterministicSteady()
 {
   JournalRecordPair pa(journal);
-  pa << "Non-linear solver for deterministic steady state By-passed " << endrec;
+  pa << "Non-linear solver for deterministic steady state skipped" << endrec;
 }
 
 void
-KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false)
+KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx)
 {
   // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals
   out.zeros();
@@ -94,16 +63,12 @@ KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx
 
 void
 KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
-                            const ConstVector &yyp, const Vector &xx) noexcept(false)
+                            const ConstVector &yyp, const Vector &xx)
 {
   // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals
   out.zeros();
 }
 
-/************************************************
- * this is main derivative calculation functin that indirectly calls dynamic.dll
- * which performs actual calculation and reorders
- ***************************************************/
 void
 KordpDynare::calcDerivativesAtSteady()
 {
@@ -136,20 +101,17 @@ KordpDynare::calcDerivativesAtSteady()
       dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p.get(), g2p.get(), g3p.get());
     }
 
-  populateDerivativesContainer(*g1p, 1, JacobianIndices);
+  populateDerivativesContainer(*g1p, 1);
 
   if (nOrder > 1)
-    populateDerivativesContainer(*g2p, 2, JacobianIndices);
+    populateDerivativesContainer(*g2p, 2);
 
   if (nOrder > 2)
-    populateDerivativesContainer(*g3p, 3, JacobianIndices);
+    populateDerivativesContainer(*g3p, 3);
 }
 
-/*******************************************************************************
- * populateDerivatives to sparse Tensor and fit it in the Derivatives Container
- *******************************************************************************/
 void
-KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector<int> &vOrder)
+KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord)
 {
   // model derivatives FSSparseTensor instance
   auto mdTi = std::make_unique<FSSparseTensor>(ord, nJcols, nY);
@@ -161,11 +123,7 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
       {
         for (int j = 0; j < g.nrows(); j++)
           {
-            double x;
+            double x = g.get(j, dynppToDyn[s[0]]);
-            if (s[0] < nJcols-nExog)
-              x = g.get(j, vOrder[s[0]]);
-            else
-              x = g.get(j, s[0]);
             if (x != 0.0)
               mdTi->insert(s, j, x);
           }
@@ -173,26 +131,16 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
       }
   else if (ord == 2)
     {
-      int nJcols1 = nJcols-nExog;
-      std::vector<int> revOrder(nJcols1);
-      for (int i = 0; i < nJcols1; i++)
-        revOrder[vOrder[i]] = i;
       for (int i = 0; i < g.nrows(); i++)
         {
           int j = static_cast<int>(g.get(i, 0))-1; // hessian indices start with 1
-          int i1 = static_cast<int>(g.get(i, 1)) -1;
+          int i1 = static_cast<int>(g.get(i, 1))-1;
           if (j < 0 || i1 < 0)
             continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
           int s0 = i1 / nJcols;
           int s1 = i1 % nJcols;
-          if (s0 < nJcols1)
+          s[0] = dynToDynpp[s0];
-            s[0] = revOrder[s0];
+          s[1] = dynToDynpp[s1];
-          else
-            s[0] = s0;
-          if (s1 < nJcols1)
-            s[1] = revOrder[s1];
-          else
-            s[1] = s1;
           if (s[1] >= s[0])
             {
               double x = g.get(i, 2);
@@ -202,67 +150,49 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const st
     }
   else if (ord == 3)
     {
-      int nJcols1 = nJcols-nExog;
       int nJcols2 = nJcols*nJcols;
-      std::vector<int> revOrder(nJcols1);
-      for (int i = 0; i < nJcols1; i++)
-        revOrder[vOrder[i]] = i;
       for (int i = 0; i < g.nrows(); i++)
         {
           int j = static_cast<int>(g.get(i, 0))-1;
-          int i1 = static_cast<int>(g.get(i, 1)) -1;
+          int i1 = static_cast<int>(g.get(i, 1))-1;
           if (j < 0 || i1 < 0)
             continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
           int s0 = i1 / nJcols2;
           int i2 = i1 % nJcols2;
           int s1 = i2 / nJcols;
           int s2 = i2 % nJcols;
-          if (s0 < nJcols1)
+          s[0] = dynToDynpp[s0];
-            s[0] = revOrder[s0];
+          s[1] = dynToDynpp[s1];
-          else
+          s[2] = dynToDynpp[s2];
-            s[0] = s0;
-          if (s1 < nJcols1)
-            s[1] = revOrder[s1];
-          else
-            s[1] = s1;
-          if (s2 < nJcols1)
-            s[2] = revOrder[s2];
-          else
-            s[2] = s2;
-          double x = g.get(i, 2);
           if (s.isSorted())
-            mdTi->insert(s, j, x);
+            {
+              double x = g.get(i, 2);
+              mdTi->insert(s, j, x);
+            }
         }
     }
 
-  // md container
-  md.remove(Symmetry{ord});
   md.insert(std::move(mdTi));
 }
 
-/*********************************************************
+/* Returns ySteady extended with leads and lags suitable for passing to
- * LLxSteady()
+   <model>_dynamic */
- * returns ySteady extended with leads and lags suitable for
- * passing to <model>_dynamic DLL
- *************************************************************/
 void
-KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false)
+KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady)
 {
-  if ((nJcols-nExog) == yS.length())
+  if (yS.length() == nJcols-nExog)
     throw DynareException(__FILE__, __LINE__, "ySteady already of right size");
 
-  // create temporary square 2D matrix size nEndo x nEndo (sparse)
+  /* Create temporary square 2D matrix size nEndo×nEndo (sparse)
-  // for the lag, current and lead blocks of the jacobian
+     for the lag, current and lead blocks of the jacobian */
   if (llxSteady.length() != nJcols-nExog)
     throw DynareException(__FILE__, __LINE__, "llxSteady has wrong size");
 
   for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
     // populate (non-sparse) vector with ysteady values
     for (int i = 0; i < nY; i++)
-      {
+      if (ll_Incidence.get(ll_row, i))
-        if (ll_Incidence.get(ll_row, i))
+        llxSteady[static_cast<int>(ll_Incidence.get(ll_row, i))-1] = yS[i];
-          llxSteady[static_cast<int>(ll_Incidence.get(ll_row, i))-1] = yS[i];
-      }
 }
 
 /************************************
@@ -291,26 +221,26 @@ KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false)
 ************************************/
 
 void
-KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
+KordpDynare::computeJacobianPermutation(const std::vector<int> &var_order)
 {
-  // create temporary square 2D matrix size nEndo x nEndo (sparse)
+  dynppToDyn.resize(nJcols);
+  // create temporary square 2D matrix size nEndo×nEndo (sparse)
   // for the lag, current and lead blocks of the jacobian
-  JacobianIndices.resize(nJcols);
   std::vector<int> tmp(nY);
-  int i, j, rjoff = nJcols-nExog-1;
+  int rjoff = nJcols-nExog-1;
 
   for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
     {
       // reorder in orde-var order & populate temporary nEndo (sparse) vector with
       // the lag, current and lead blocks of the jacobian respectively
-      for (i = 0; i < nY; i++)
+      for (int i = 0; i < nY; i++)
-        tmp[i] = static_cast<int>(ll_Incidence.get(ll_row, varOrder[i]-1));
+        tmp[i] = static_cast<int>(ll_Incidence.get(ll_row, var_order[i]-1));
       // write the reordered blocks back to the jacobian
       // in reverse order
-      for (j = nY-1; j >= 0; j--)
+      for (int j = nY-1; j >= 0; j--)
         if (tmp[j])
           {
-            JacobianIndices[rjoff] = tmp[j] -1;
+            dynppToDyn[rjoff] = tmp[j]-1;
             rjoff--;
             if (rjoff < 0)
               break;
@@ -318,13 +248,15 @@ KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
     }
 
   //add the indices for the nExog exogenous jacobians
-  for (j = nJcols-nExog; j < nJcols; j++)
+  for (int j = nJcols-nExog; j < nJcols; j++)
-    JacobianIndices[j] = j;
+    dynppToDyn[j] = j;
+
+  // Compute reverse ordering
+  dynToDynpp.resize(nJcols);
+  for (int i = 0; i < nJcols; i++)
+    dynToDynpp[dynppToDyn[i]] = i;
 }
 
-/**************************************************************************************/
-/*       DynareNameList class                                                         */
-/**************************************************************************************/
 
 DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg)
   : names(std::move(names_arg))
@@ -334,7 +266,7 @@ DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::strin
 DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const DynareNameList &dnl,
                                          const DynareNameList &denl)
 {
-  for (int i = 0; i < dynare.nys(); i++)
+  for (int i = 0; i < dynare.nYs; i++)
     names.emplace_back(dnl.getName(i+dynare.nstat()));
   for (int i = 0; i < dynare.nexog(); i++)
     names.emplace_back(denl.getName(i));

mex/sources/k_order_perturbation/k_ord_dynare.hh

@@ -21,39 +21,33 @@
 #define K_ORD_DYNARE3_H
 
 #include <vector>
+#include <string>
 #include <memory>
 
+#include "journal.hh"
+#include "Vector.hh"
+#include "twod_matrix.hh"
 #include "t_container.hh"
 #include "sparse_tensor.hh"
 #include "decision_rule.hh"
 #include "dynamic_model.hh"
-
+#include "dynamic_abstract_class.hh"
-#include "exception.hh"
-#include "dynare_exception.hh"
-#include "fs_tensor.hh"
-#include "SylvException.hh"
-#include "tl_exception.hh"
-#include "kord_exception.hh"
-#include "nlsolve.hh"
-#include "approximation.hh"
 
 class KordpDynare;
 
-/*////////////////////////////////////////////*/
+// instantiations of pure abstract class NameList in dynamic_model.h
-// instantiations of pure abstract class NameList in dynamic_model.h:
-/*////////////////////////////////////////////*/
 class DynareNameList : public NameList
 {
   std::vector<std::string> names;
 public:
   DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg);
   int
-  getNum() const
+  getNum() const override
   {
     return static_cast<int>(names.size());
   }
   const std::string &
-  getName(int i) const
+  getName(int i) const override
   {
     return names[i];
   }
@@ -66,44 +60,37 @@ public:
   DynareStateNameList(const KordpDynare &dynare, const DynareNameList &dnl,
                       const DynareNameList &denl);
   int
-  getNum() const
+  getNum() const override
   {
     return static_cast<int>(names.size());
   }
   const std::string &
-  getName(int i) const
+  getName(int i) const override
   {
     return names[i];
   }
 };
-/*********************************************/
+
 // The following only implements DynamicModel with help of ogdyn::DynareModel
 // instantiation of pure abstract DynamicModel decl. in dynamic_model.h
-class DynamicModelAC;
-class DynamicModelDLL;
-class DynamicModelMFile;
-
 class KordpDynare : public DynamicModel
 {
-  friend class DynareNameList;
+public:
-  friend class DynareStateNameList;
-  friend class DynamicModelDLL;
-  friend class DynamicModelMFile;
-
   const int nStat;
   const int nBoth;
   const int nPred;
   const int nForw;
   const int nExog;
   const int nPar;
-  const int nYs; // ={npred + nboth ; }
+  const int nYs;      // = npred + nboth
-  const int nYss; // nyss ={ nboth + nforw ; }
+  const int nYss;     // = nboth + nforw
-  const int nY;  // = num_endo={ nstat + npred + nboth + nforw ; }
+  const int nY;       // = nstat + npred + nboth + nforw
-  const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw
+  const int nJcols;   // nb of jacobian columns = nExog+nY+nYs+nYss
-  const Vector &NNZD;  /* the total number of non-zero derivative elements
+  const Vector &NNZD; /* the total number of non-zero derivative elements
-                          where hessian is 2nd : NZZD(order=2) */
+                         where hessian is 2nd : NNZD(order=2) */
   const int nSteps;
   const int nOrder;
+private:
   Journal &journal;
   Vector &ySteady;
   Vector &params;
@@ -111,141 +98,110 @@ class KordpDynare : public DynamicModel
   TensorContainer<FSSparseTensor> md; // ModelDerivatives
   DynareNameList dnl, denl;
   DynareStateNameList dsnl;
-  const double ss_tol;
-  const std::vector<int> &varOrder;
   const TwoDMatrix &ll_Incidence;
-  double qz_criterium;
+  std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
-  std::vector<int> JacobianIndices;
+  std::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
 
   std::unique_ptr<TwoDMatrix> g1p, g2p, g3p;
+
+  std::unique_ptr<DynamicModelAC> dynamicModelFile;
 public:
-  KordpDynare(const std::vector<std::string> &endo, int num_endo,
+  KordpDynare(const std::vector<std::string> &endo,
               const std::vector<std::string> &exo, int num_exo, int num_par,
               Vector &ySteady, TwoDMatrix &vCov, Vector &params, int nstat, int nPred,
-              int nforw, int nboth, const int nJcols, const Vector &NNZD,
+              int nforw, int nboth, const Vector &NNZD,
-              const int nSteps, const int ord,
+              int nSteps, int ord,
-              Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
+              Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
               const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
-              double qz_criterium) noexcept(false);
+              std::unique_ptr<TwoDMatrix> g1_arg,
-  KordpDynare(const std::vector<std::string> &endo, int num_endo,
+              std::unique_ptr<TwoDMatrix> g2_arg, std::unique_ptr<TwoDMatrix> g3_arg);
-              const std::vector<std::string> &exo, int num_exo, int num_par,
-              Vector &ySteady, TwoDMatrix &vCov, Vector &params, int nstat, int nPred,
-              int nforw, int nboth, const int nJcols, const Vector &NNZD,
-              const int nSteps, const int ord,
-              Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg, double sstol,
-              const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
-              double qz_criterium, std::unique_ptr<TwoDMatrix> g1_arg,
-              std::unique_ptr<TwoDMatrix> g2_arg, std::unique_ptr<TwoDMatrix> g3_arg) noexcept(false);
 
   int
-  nstat() const
+  nstat() const override
   {
     return nStat;
   }
   int
-  nboth() const
+  nboth() const override
   {
     return nBoth;
   }
   int
-  npred() const
+  npred() const override
   {
     return nPred;
   }
   int
-  nforw() const
+  nforw() const override
   {
     return nForw;
   }
   int
-  nexog() const
+  nexog() const override
   {
     return nExog;
   }
   int
-  nys() const
+  order() const override
-  {
-    return nYs;
-  }
-  int
-  nyss() const
-  {
-    return nYss;
-  }
-  int
-  ny() const
-  {
-    return nY;
-  }
-  int
-  steps() const
-  {
-    return nSteps;
-  }
-  int
-  order() const
   {
     return nOrder;
   }
   const NameList &
-  getAllEndoNames() const
+  getAllEndoNames() const override
   {
     return dnl;
   }
   const NameList &
-  getStateNames() const
+  getStateNames() const override
   {
     return dsnl;
   }
   const NameList &
-  getExogNames() const
+  getExogNames() const override
   {
     return denl;
   }
   const TwoDMatrix &
-  getVcov() const
+  getVcov() const override
   {
     return vCov;
   }
-  Vector &
-  getParams()
-  {
-    return params;
-  }
 
   const TensorContainer<FSSparseTensor> &
-  getModelDerivatives() const
+  getModelDerivatives() const override
   {
     return md;
   }
   const Vector &
-  getSteady() const
+  getSteady() const override
   {
     return ySteady;
   }
   Vector &
-  getSteady()
+  getSteady() override
   {
     return ySteady;
   }
 
-  void solveDeterministicSteady();
+  void solveDeterministicSteady() override;
-  void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false);
+  void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) override;
   void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
-                      const ConstVector &yyp, const Vector &xx) noexcept(false);
+                      const ConstVector &yyp, const Vector &xx) override;
-  void calcDerivativesAtSteady();
+  void calcDerivativesAtSteady() override;
-  std::unique_ptr<DynamicModelAC> dynamicModelFile;
   std::unique_ptr<DynamicModel>
-  clone() const
+  clone() const override
   {
     std::cerr << "KordpDynare::clone() not implemented" << std::endl;
     exit(EXIT_FAILURE);
   }
-  void LLxSteady(const Vector &yS, Vector &llxSteady) noexcept(false); // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags
-
 private:
-  void ReorderDynareJacobianIndices() noexcept(false);
+  // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags
-  void populateDerivativesContainer(const TwoDMatrix &g, int ord, const std::vector<int> &vOrder);
+  void LLxSteady(const Vector &yS, Vector &llxSteady);
+  /* Computes the permutations mapping back and forth between Dynare and
+     Dynare++ orderings of variables */
+  void computeJacobianPermutation(const std::vector<int> &var_order);
+  // Fills tensor container (with reordering) at a given order
+  void populateDerivativesContainer(const TwoDMatrix &g, int ord);
 };
 
 #endif

mex/sources/k_order_perturbation/k_order_perturbation.cc

@@ -17,36 +17,33 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-/*
+/* Defines the entry point for the k-order perturbation application DLL.
-  Defines the entry point for the k-order perturbation application DLL.
 
-  Inputs:
+   See matlab/mex/k_order_perturbation.m for a description of inputs and
-  1) dr
+   outputs.
-  2) M_
-  3) options
-
-  Outputs:
-  - if order == 1: only g_1
-  - if order == 2: g_0, g_1, g_2
-  - if order == 3: g_0, g_1, g_2, g_3
 */
 
 #include "dynamic_m.hh"
 #include "dynamic_dll.hh"
+#include "k_ord_dynare.hh"
+
+#include "approximation.hh"
+#include "exception.hh"
+#include "dynare_exception.hh"
+#include "kord_exception.hh"
+#include "tl_exception.hh"
+#include "SylvException.hh"
 
 #include <algorithm>
-#include <cmath>
 #include <cassert>
 
 #include "dynmex.h"
 
-//////////////////////////////////////////////////////
+/* Convert MATLAB Dynare endo and exo names array to a vector<string> array of
-// Convert MATLAB Dynare endo and exo names array to a vector<string> array of string pointers
+   string pointers. MATLAB “mx” function returns a long string concatenated by
-// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
+   columns rather than rows hence a rather low level approach is needed. */
-// hence a rather low level approach is needed
-///////////////////////////////////////////////////////
 void
-DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std::vector<std::string> &out)
+DynareMxArrayToString(const mxArray *mxFldp, int len, int width, std::vector<std::string> &out)
 {
   char *cNamesCharStr = mxArrayToString(mxFldp);
 
@@ -55,7 +52,7 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std
   for (int i = 0; i < width; i++)
     for (int j = 0; j < len; j++)
       // Allow alphanumeric and underscores "_" only:
-      if (isalnum(cNamesCharStr[j+i*len]) || (cNamesCharStr[j+i*len] == '_'))
+      if (std::isalnum(cNamesCharStr[j+i*len]) || (cNamesCharStr[j+i*len] == '_'))
         out[j] += cNamesCharStr[j+i*len];
 }
 
@@ -122,14 +119,10 @@ extern "C" {
     const int nStat = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "npred");
     const int nPred = static_cast<int>(mxGetScalar(mxFldp));
-    mxFldp = mxGetField(M_, 0, "nspred");
-    const int nsPred = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nboth");
     const int nBoth = static_cast<int>(mxGetScalar(mxFldp));
     mxFldp = mxGetField(M_, 0, "nfwrd");
     const int nForw = static_cast<int>(mxGetScalar(mxFldp));
-    mxFldp = mxGetField(M_, 0, "nsfwrd");
-    const int nsForw = static_cast<int>(mxGetScalar(mxFldp));
 
     mxFldp = mxGetField(M_, 0, "exo_nbr");
     const int nExog = static_cast<int>(mxGetScalar(mxFldp));
@@ -158,13 +151,10 @@ extern "C" {
       DYN_MEX_FUNC_ERR_MSG_TXT(("dynare:k_order_perturbation: Incorrect length of lead lag incidences: ncol="
                                 + std::to_string(npar) + " != nEndo=" + std::to_string(nEndo)).c_str());
 
-    //get NNZH =NNZD(2) = the total number of non-zero Hessian elements
     mxFldp = mxGetField(M_, 0, "NNZDerivatives");
     Vector NNZD{mxFldp};
     if (NNZD[kOrder-1] == -1)
-      DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command");
+      DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the 'stoch_simul' command");
-
-    const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
 
     mxFldp = mxGetField(M_, 0, "var_order_endo_names");
     const int nendo = static_cast<int>(mxGetM(mxFldp));
@@ -182,9 +172,9 @@ extern "C" {
       DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input parameters.");
 
     std::unique_ptr<TwoDMatrix> g1m, g2m, g3m;
-    // derivatives passed as arguments */
     if (nrhs > 3)
       {
+        // Derivatives have been passed as arguments
         const mxArray *g1 = prhs[3];
         int m = static_cast<int>(mxGetM(g1));
         int n = static_cast<int>(mxGetN(g1));
@@ -206,7 +196,6 @@ extern "C" {
       }
 
     const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
-    const double sstol = 1.e-13; //NL solver tolerance from
 
     try
       {
@@ -222,10 +211,10 @@ extern "C" {
         TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);
 
         // make KordpDynare object
-        KordpDynare dynare(endoNames, nEndo, exoNames, nExog, nPar,
+        KordpDynare dynare(endoNames, exoNames, nExog, nPar,
                            ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
-                           jcols, NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
+                           NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
-                           sstol, var_order_vp, llincidence, qz_criterium,
+                           var_order_vp, llincidence,
                            std::move(g1m), std::move(g2m), std::move(g3m));
 
         // construct main K-order approximation class