diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 7a4436eff95bf7d840da3b36949bf94445eb98e2..90ee33610b0a75fa3f033e934092b31a481eb2ad 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2008-2022 Dynare Team
+ * Copyright © 2008-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -161,15 +161,6 @@ extern "C"
const int nEndo = get_int_field(M_mx, "endo_nbr");
const int nPar = get_int_field(M_mx, "param_nbr");
- const mxArray* lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
- if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx)
- && !mxIsComplex(lead_lag_incidence_mx) && !mxIsSparse(lead_lag_incidence_mx)
- && mxGetM(lead_lag_incidence_mx) == 3
- && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(nEndo)))
- mexErrMsgTxt("M_.lead_lag_incidence should be a real dense matrix with 3 rows and "
- "M_.endo_nbr columns");
- ConstTwoDMatrix llincidence {lead_lag_incidence_mx};
-
const mxArray* nnzderivatives_mx = mxGetField(M_mx, 0, "NNZDerivatives");
if (!(nnzderivatives_mx && mxIsDouble(nnzderivatives_mx) && !mxIsComplex(nnzderivatives_mx)
&& !mxIsSparse(nnzderivatives_mx)))
@@ -217,6 +208,31 @@ extern "C"
std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx) + nEndo, dr_order.begin(),
[](double x) { return static_cast<int>(x) - 1; });
+ // Retrieve sparse indices for dynamic model
+
+ const mxArray* dynamic_g1_sparse_rowval_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_rowval");
+ if (!(dynamic_g1_sparse_rowval_mx && mxIsInt32(dynamic_g1_sparse_rowval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_rowval should be an int32 array");
+
+ const mxArray* dynamic_g1_sparse_colval_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_colval");
+ if (!(dynamic_g1_sparse_colval_mx && mxIsInt32(dynamic_g1_sparse_colval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colval should be an int32 array");
+
+ const mxArray* dynamic_g1_sparse_colptr_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_colptr");
+ if (!(dynamic_g1_sparse_colptr_mx && mxIsInt32(dynamic_g1_sparse_colptr_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colptr should be an int32 array");
+
+ std::vector<const mxArray*> dynamic_gN_sparse_indices;
+ for (int o {2}; o <= kOrder; o++)
+ {
+ using namespace std::string_literals;
+ auto fieldname {"dynamic_g"s + std::to_string(o) + "_sparse_indices"};
+ const mxArray* indices = mxGetField(M_mx, 0, fieldname.c_str());
+ if (!(indices && mxIsInt32(indices)))
+ mexErrMsgTxt(("M_."s + fieldname + " should be an int32 array").c_str());
+ dynamic_gN_sparse_indices.push_back(indices);
+ }
+
const int nSteps
= 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
@@ -229,9 +245,13 @@ extern "C"
std::unique_ptr<DynamicModelAC> dynamicModelFile;
if (use_dll)
- dynamicModelFile = std::make_unique<DynamicModelDLL>(fname, ntt, kOrder);
+ dynamicModelFile = std::make_unique<DynamicModelDLL>(
+ fname, kOrder, dynamic_g1_sparse_rowval_mx, dynamic_g1_sparse_colval_mx,
+ dynamic_g1_sparse_colptr_mx, dynamic_gN_sparse_indices, ntt);
else
- dynamicModelFile = std::make_unique<DynamicModelMFile>(fname, ntt);
+ dynamicModelFile = std::make_unique<DynamicModelMFile>(
+ fname, kOrder, dynamic_g1_sparse_rowval_mx, dynamic_g1_sparse_colval_mx,
+ dynamic_g1_sparse_colptr_mx, dynamic_gN_sparse_indices);
// intiate tensor library
TLStatic::init(kOrder, nStat + 2 * nPred + 3 * nBoth + 2 * nForw + nExog);
@@ -242,7 +262,7 @@ extern "C"
// make KordpDynare object
KordpDynare dynare(endoNames, exoNames, nExog, nPar, ySteady, vCov, modParams, nStat, nPred,
nForw, nBoth, NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
- dr_order, llincidence);
+ dr_order);
// construct main K-order approximation class
Approximation app(dynare, journal, nSteps, false, pruning, qz_criterium);
diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc
index 130afa58a9e9f4314a187056c999b38471478989..9f0c72a6f3d3be2e03c8e666404b394271d289ba 100644
--- a/mex/sources/k_order_welfare/k_order_welfare.cc
+++ b/mex/sources/k_order_welfare/k_order_welfare.cc
@@ -199,14 +199,6 @@ extern "C"
const int nEndo = get_int_field(M_mx, "endo_nbr");
const int nPar = get_int_field(M_mx, "param_nbr");
- const mxArray* lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
- if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx)
- && mxGetM(lead_lag_incidence_mx) == 3
- && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(nEndo)))
- mexErrMsgTxt("M_.lead_lag_incidence should be a double precision matrix with 3 rows and "
- "M_.endo_nbr columns");
- ConstTwoDMatrix llincidence {lead_lag_incidence_mx};
-
const mxArray* nnzderivatives_mx = mxGetField(M_mx, 0, "NNZDerivatives");
if (!(nnzderivatives_mx && mxIsDouble(nnzderivatives_mx)))
mexErrMsgTxt("M_.NNZDerivatives should be a double precision array");
@@ -253,6 +245,33 @@ extern "C"
std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx) + nEndo, dr_order.begin(),
[](double x) { return static_cast<int>(x) - 1; });
+ // Retrieve sparse indices for dynamic model
+
+ const mxArray* dynamic_g1_sparse_rowval_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_rowval");
+ if (!(dynamic_g1_sparse_rowval_mx && mxIsInt32(dynamic_g1_sparse_rowval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_rowval should be an int32 array");
+
+ const mxArray* dynamic_g1_sparse_colval_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_colval");
+ if (!(dynamic_g1_sparse_colval_mx && mxIsInt32(dynamic_g1_sparse_colval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colval should be an int32 array");
+
+ const mxArray* dynamic_g1_sparse_colptr_mx = mxGetField(M_mx, 0, "dynamic_g1_sparse_colptr");
+ if (!(dynamic_g1_sparse_colptr_mx && mxIsInt32(dynamic_g1_sparse_colptr_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colptr should be an int32 array");
+
+ std::vector<const mxArray*> dynamic_gN_sparse_indices;
+ for (int o {2}; o <= kOrder; o++)
+ {
+ using namespace std::string_literals;
+ auto fieldname {"dynamic_g"s + std::to_string(o) + "_sparse_indices"};
+ const mxArray* indices = mxGetField(M_mx, 0, fieldname.c_str());
+ if (!(indices && mxIsInt32(indices)))
+ mexErrMsgTxt(("M_."s + fieldname + " should be an int32 array").c_str());
+ dynamic_gN_sparse_indices.push_back(indices);
+ }
+
+ // Retrieve sparse indices for objective model
+
const mxArray* objective_g1_sparse_rowval_mx
= mxGetField(M_mx, 0, "objective_g1_sparse_rowval");
if (!(objective_g1_sparse_rowval_mx && mxIsInt32(objective_g1_sparse_rowval_mx)))
@@ -289,9 +308,13 @@ extern "C"
std::unique_ptr<DynamicModelAC> dynamicModelFile;
if (use_dll)
- dynamicModelFile = std::make_unique<DynamicModelDLL>(fname, ntt, kOrder);
+ dynamicModelFile = std::make_unique<DynamicModelDLL>(
+ fname, kOrder, dynamic_g1_sparse_rowval_mx, dynamic_g1_sparse_colval_mx,
+ dynamic_g1_sparse_colptr_mx, dynamic_gN_sparse_indices, ntt);
else
- dynamicModelFile = std::make_unique<DynamicModelMFile>(fname, ntt);
+ dynamicModelFile = std::make_unique<DynamicModelMFile>(
+ fname, kOrder, dynamic_g1_sparse_rowval_mx, dynamic_g1_sparse_colval_mx,
+ dynamic_g1_sparse_colptr_mx, dynamic_gN_sparse_indices);
// intiate tensor library
TLStatic::init(kOrder, nStat + 2 * nPred + 3 * nBoth + 2 * nForw + nExog);
@@ -302,7 +325,7 @@ extern "C"
// make KordpDynare object
KordpDynare dynare(endoNames, exoNames, nExog, nPar, ySteady, vCov, modParams, nStat, nPred,
nForw, nBoth, NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
- dr_order, llincidence);
+ dr_order);
// construct main K-order approximation class
Approximation app(dynare, journal, nSteps, false, pruning, qz_criterium);
diff --git a/mex/sources/libkorder/dynamic_abstract_class.hh b/mex/sources/libkorder/dynamic_abstract_class.hh
index 951e845bfff3284255f9e704b3810e7ce49bd976..1a1aa9a175b0313340899c3b24a9d710881b8716 100644
--- a/mex/sources/libkorder/dynamic_abstract_class.hh
+++ b/mex/sources/libkorder/dynamic_abstract_class.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -22,17 +22,36 @@
#include <vector>
-#include "twod_matrix.hh"
+#include "dynmex.h"
+
+#include "Vector.hh"
+#include "sparse_tensor.hh"
+#include "t_container.hh"
class DynamicModelAC
{
protected:
- int ntt; // Size of vector of temporary terms
+ const int order;
+ const mxArray *const dynamic_g1_sparse_rowval_mx, *const dynamic_g1_sparse_colval_mx,
+ *const dynamic_g1_sparse_colptr_mx;
+ // Stores M_.dynamic_gN_sparse_indices, starting from N=2
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices;
+
public:
- DynamicModelAC(int ntt_arg) : ntt {ntt_arg} {};
+ DynamicModelAC(int order_arg, const mxArray* dynamic_g1_sparse_rowval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colptr_mx_arg,
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices_arg) :
+ order {order_arg},
+ dynamic_g1_sparse_rowval_mx {dynamic_g1_sparse_rowval_mx_arg},
+ dynamic_g1_sparse_colval_mx {dynamic_g1_sparse_colval_mx_arg},
+ dynamic_g1_sparse_colptr_mx {dynamic_g1_sparse_colptr_mx_arg},
+ dynamic_gN_sparse_indices {dynamic_gN_sparse_indices_arg} {};
virtual ~DynamicModelAC() = default;
virtual void eval(const Vector& y, const Vector& x, const Vector& params, const Vector& ySteady,
- Vector& residual, std::vector<TwoDMatrix>& md)
+ Vector& residual, const std::map<int, int>& dynToDynpp,
+ TensorContainer<FSSparseTensor>& derivatives) noexcept(false)
= 0;
};
+
#endif
diff --git a/mex/sources/libkorder/dynamic_dll.cc b/mex/sources/libkorder/dynamic_dll.cc
index 960f63a48e9d124dd8a38a9f2a6ef5cacba20de2..32d1cbe7855effc13c504145081f0a873a67cc8a 100644
--- a/mex/sources/libkorder/dynamic_dll.cc
+++ b/mex/sources/libkorder/dynamic_dll.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2008-2023 Dynare Team
+ * Copyright © 2008-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -22,76 +22,159 @@
#include <cassert>
#include <iostream>
-DynamicModelDLL::DynamicModelDLL(const std::string& modName, int ntt_arg, int order) :
- DynamicModelAC(ntt_arg)
+DynamicModelDLL::DynamicModelDLL(const std::string& modName, int order_arg,
+ const mxArray* dynamic_g1_sparse_rowval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colptr_mx_arg,
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices_arg,
+ int ntt) :
+ DynamicModelAC {order_arg, dynamic_g1_sparse_rowval_mx_arg, dynamic_g1_sparse_colval_mx_arg,
+ dynamic_g1_sparse_colptr_mx_arg, dynamic_gN_sparse_indices_arg},
+ dynamic_fcts(order + 1),
+ dynamic_fcts_tt(order + 1),
+ tt_tmp(ntt),
+ derivatives_tmp(order)
{
- std::string fName;
+ using namespace std::string_literals;
+
+ for (int o {0}; o <= order; o++)
+ {
+ std::string func_name {"dynamic_"s + (o == 0 ? "resid" : "g"s + std::to_string(o))};
+ std::string mex_filename
+ {
#if !defined(__CYGWIN32__) && !defined(_WIN32)
- fName = "./";
+ "./"s +
#endif
- fName += "+" + modName + "/dynamic" + MEXEXT;
+ "+"s + modName + "/+sparse/" + func_name + MEXEXT
+ };
-#if defined(__CYGWIN32__) || defined(_WIN32)
- dynamicHinstance = LoadLibrary(fName.c_str());
-#else // GNU/Linux or Mac
- dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW);
+ mex_handle_t handle {load_mex(mex_filename)};
+ if (handle)
+ mex_handles.push_back(handle);
+ else
+ {
+ std::ranges::for_each(mex_handles, unload_mex);
+
+ throw DynareException(__FILE__, __LINE__,
+ "Error when loading " + mex_filename
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ + ": " + dlerror()
#endif
- if (!dynamicHinstance)
- throw DynareException(__FILE__, __LINE__,
- "Error when loading " + fName
+ );
+ }
+
+ std::tie(dynamic_fcts[o], dynamic_fcts_tt[o])
+ = getSymbolsFromDLL(func_name, mex_handles.back());
+ if (!dynamic_fcts[o] || !dynamic_fcts_tt[o])
+ {
+ std::ranges::for_each(mex_handles, unload_mex);
+
+ throw DynareException(__FILE__, __LINE__,
+ "Error when loading symbols from " + mex_filename
#if !defined(__CYGWIN32__) && !defined(_WIN32)
- + ": " + dlerror()
+ + ": " + dlerror()
#endif
- );
+ );
+ }
+ }
- dynamic_tt.resize(order + 1);
+ derivatives_tmp[0].resize(mxGetNumberOfElements(dynamic_g1_sparse_rowval_mx));
+ for (int o {2}; o <= order; o++)
+ derivatives_tmp[o - 1].resize(mxGetM(dynamic_gN_sparse_indices[o - 2]));
+}
- std::tie(dynamic_resid, dynamic_tt[0])
- = getSymbolsFromDLL<dynamic_resid_or_g1_fct>("dynamic_resid", fName);
- std::tie(dynamic_g1, dynamic_tt[1])
- = getSymbolsFromDLL<dynamic_resid_or_g1_fct>("dynamic_g1", fName);
+DynamicModelDLL::~DynamicModelDLL()
+{
+ std::ranges::for_each(mex_handles, unload_mex);
+}
- dynamic_higher_deriv.resize(std::max(0, order - 1));
- for (int i = 2; i <= order; i++)
- std::tie(dynamic_higher_deriv[i - 2], dynamic_tt[i])
- = getSymbolsFromDLL<dynamic_higher_deriv_fct>("dynamic_g" + std::to_string(i), fName);
+DynamicModelDLL::mex_handle_t
+DynamicModelDLL::load_mex(const std::string& mex_filename)
+{
+#if defined(__CYGWIN32__) || defined(_WIN32)
+ return LoadLibrary(mex_filename.c_str());
+#else // GNU/Linux or Mac
+ return dlopen(mex_filename.c_str(), RTLD_NOW);
+#endif
+}
- tt.resize(ntt);
+void
+DynamicModelDLL::unload_mex(mex_handle_t handle)
+{
+#if defined(__CYGWIN32__) || defined(_WIN32)
+ FreeLibrary(handle);
+#else
+ dlclose(handle);
+#endif
}
-DynamicModelDLL::~DynamicModelDLL()
+std::pair<dynamic_fct, dynamic_tt_fct>
+DynamicModelDLL::getSymbolsFromDLL(const std::string& func_name, mex_handle_t handle)
{
#if defined(__CYGWIN32__) || defined(_WIN32)
- auto result = FreeLibrary(dynamicHinstance);
- if (result == 0)
- {
- std::cerr << "Can't free the *_dynamic DLL" << std::endl;
- exit(EXIT_FAILURE);
- }
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Wcast-function-type"
+ return {reinterpret_cast<dynamic_fct>(GetProcAddress(handle, func_name.c_str())),
+ reinterpret_cast<dynamic_tt_fct>(GetProcAddress(handle, (func_name + "_tt").c_str()))};
+# pragma GCC diagnostic pop
#else
- dlclose(dynamicHinstance);
+ return {reinterpret_cast<dynamic_fct>(dlsym(handle, func_name.c_str())),
+ reinterpret_cast<dynamic_tt_fct>(dlsym(handle, (func_name + "_tt").c_str()))};
#endif
}
void
DynamicModelDLL::eval(const Vector& y, const Vector& x, const Vector& modParams,
- const Vector& ySteady, Vector& residual,
- std::vector<TwoDMatrix>& md) noexcept(false)
+ const Vector& ySteady, Vector& residual, const std::map<int, int>& dynToDynpp,
+ TensorContainer<FSSparseTensor>& derivatives) noexcept(false)
{
- assert(md.size() == dynamic_tt.size() - 1);
-
- for (size_t i = 0; i < dynamic_tt.size(); i++)
+ for (int o {0}; o <= order; o++)
{
- dynamic_tt[i](y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.data());
- if (i == 0)
- dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.data(),
- residual.base());
- else if (i == 1)
- dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.data(),
- md[0].base());
- else
- dynamic_higher_deriv[i - 2](y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0,
- tt.data(), &md[i - 1].get(0, 0), &md[i - 1].get(0, 1),
- &md[i - 1].get(0, 2));
+ dynamic_fcts_tt[o](y.base(), x.base(), modParams.base(), ySteady.base(), tt_tmp.data());
+ dynamic_fcts[o](y.base(), x.base(), modParams.base(), ySteady.base(), tt_tmp.data(),
+ o == 0 ? residual.base() : derivatives_tmp[o - 1].data());
+
+ if (o >= 1)
+ {
+ IntSequence s(o, 0);
+ auto tensor = std::make_unique<FSSparseTensor>(o, dynToDynpp.size(), residual.length());
+ size_t nnz {derivatives_tmp[o - 1].size()};
+
+ for (size_t k {0}; k < nnz; k++)
+ {
+ int rowval;
+ if (o == 1)
+ {
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T* sparse_rowval {mxGetInt32s(dynamic_g1_sparse_rowval_mx)};
+ const int32_T* sparse_colval {mxGetInt32s(dynamic_g1_sparse_colval_mx)};
+#else
+ const int32_T* sparse_rowval {
+ static_cast<const int32_T*>(mxGetData(dynamic_g1_sparse_rowval_mx))};
+ const int32_T* sparse_colval {
+ static_cast<const int32_T*>(mxGetData(dynamic_g1_sparse_colval_mx))};
+#endif
+ s[0] = dynToDynpp.at(sparse_colval[k] - 1);
+ rowval = sparse_rowval[k] - 1;
+ }
+ else
+ {
+ const mxArray* sparse_indices_mx {dynamic_gN_sparse_indices[o - 2]};
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T* sparse_indices {mxGetInt32s(sparse_indices_mx)};
+#else
+ const int32_T* sparse_indices {
+ static_cast<const int32_T*>(mxGetData(sparse_indices_mx))};
+#endif
+ for (int i {0}; i < o; i++)
+ s[i] = dynToDynpp.at(sparse_indices[k + (i + 1) * nnz] - 1);
+ s.sort();
+ rowval = sparse_indices[k] - 1;
+ }
+ tensor->insert(s, rowval, derivatives_tmp[o - 1][k]);
+ }
+
+ derivatives.insert(std::move(tensor));
+ }
}
}
diff --git a/mex/sources/libkorder/dynamic_dll.hh b/mex/sources/libkorder/dynamic_dll.hh
index 8eca28d9caa257a752e2694e67b663c96602eb06..9f3f365171ff5e6d420961a0ccb4ea23427ac2b9 100644
--- a/mex/sources/libkorder/dynamic_dll.hh
+++ b/mex/sources/libkorder/dynamic_dll.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2008-2023 Dynare Team
+ * Copyright © 2008-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -37,15 +37,10 @@
#include "dynamic_abstract_class.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_or_g1_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* resid_or_g1);
-using dynamic_higher_deriv_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* g_i, double* g_j, double* g_v);
+using dynamic_tt_fct = void (*)(const double* y, const double* x, const double* params,
+ const double* steady_state, double* T);
+using dynamic_fct = void (*)(const double* y, const double* x, const double* params,
+ const double* steady_state, const double* T, double* values);
/**
* creates pointer to Dynamic function inside <model>_dynamic.dll
@@ -54,56 +49,38 @@ using dynamic_higher_deriv_fct = void (*)(const double* y, const double* x, int
class DynamicModelDLL : public DynamicModelAC
{
private:
- std::vector<dynamic_tt_fct> dynamic_tt;
- dynamic_resid_or_g1_fct dynamic_resid, dynamic_g1;
- std::vector<dynamic_higher_deriv_fct> dynamic_higher_deriv; // Index 0 is g2
+ using mex_handle_t =
#if defined(_WIN32) || defined(__CYGWIN32__)
- HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
+ HINSTANCE
#else
- void* dynamicHinstance; // and in Linux or Mac
+ void*
#endif
- std::vector<double> tt; // Vector of temporary terms
-
- template<typename T>
- std::pair<T, dynamic_tt_fct>
- getSymbolsFromDLL(const std::string& funcname, const std::string& fName)
- {
- dynamic_tt_fct tt;
- T deriv;
-#if defined(__CYGWIN32__) || defined(_WIN32)
-# pragma GCC diagnostic push
-# pragma GCC diagnostic ignored "-Wcast-function-type"
- deriv = reinterpret_cast<T>(GetProcAddress(dynamicHinstance, funcname.c_str()));
- tt = reinterpret_cast<dynamic_tt_fct>(
- GetProcAddress(dynamicHinstance, (funcname + "_tt").c_str()));
-# pragma GCC diagnostic pop
-#else
- deriv = reinterpret_cast<T>(dlsym(dynamicHinstance, funcname.c_str()));
- tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, (funcname + "_tt").c_str()));
-#endif
- if (!deriv || !tt)
- {
-#if defined(__CYGWIN32__) || defined(_WIN32)
- FreeLibrary(dynamicHinstance);
-#else
- dlclose(dynamicHinstance);
-#endif
- throw DynareException(__FILE__, __LINE__,
- "Error when loading symbols from " + fName
-#if !defined(__CYGWIN32__) && !defined(_WIN32)
- + ": " + dlerror()
-#endif
- );
- }
- return {deriv, tt};
- }
+ ;
+ std::vector<mex_handle_t> mex_handles;
+ std::vector<dynamic_fct> dynamic_fcts; // Index 0 is resid
+ std::vector<dynamic_tt_fct> dynamic_fcts_tt;
+ std::vector<double> tt_tmp; // Vector of temporary terms
+ std::vector<std::vector<double>> derivatives_tmp;
public:
// construct and load Dynamic model DLL
- explicit DynamicModelDLL(const std::string& fname, int ntt_arg, int order);
+ DynamicModelDLL(const std::string& fname, int order_arg,
+ const mxArray* dynamic_g1_sparse_rowval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colptr_mx_arg,
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices_arg, int ntt);
~DynamicModelDLL() override;
void eval(const Vector& y, const Vector& x, const Vector& params, const Vector& ySteady,
- Vector& residual, std::vector<TwoDMatrix>& md) override;
+ Vector& residual, const std::map<int, int>& dynToDynpp,
+ TensorContainer<FSSparseTensor>& derivatives) noexcept(false) override;
+
+private:
+ // Returns a non-null pointer on success
+ static mex_handle_t load_mex(const std::string& mex_filename);
+ static void unload_mex(mex_handle_t handle);
+ // Returns non-null pointers on success
+ static std::pair<dynamic_fct, dynamic_tt_fct> getSymbolsFromDLL(const std::string& func_name,
+ mex_handle_t handle);
};
#endif
diff --git a/mex/sources/libkorder/dynamic_m.cc b/mex/sources/libkorder/dynamic_m.cc
index 6ad53131699a33889b88da50e93797e134aa167c..087186c5a347a6ecf8cb8489f6dc3ba40bba6cfd 100644
--- a/mex/sources/libkorder/dynamic_m.cc
+++ b/mex/sources/libkorder/dynamic_m.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -26,8 +26,13 @@
#include "dynamic_m.hh"
-DynamicModelMFile::DynamicModelMFile(const std::string& modName, int ntt_arg) :
- DynamicModelAC(ntt_arg), DynamicMFilename {modName + ".dynamic"}
+DynamicModelMFile::DynamicModelMFile(
+ const std::string& modName, int order_arg, const mxArray* dynamic_g1_sparse_rowval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colval_mx_arg, const mxArray* dynamic_g1_sparse_colptr_mx_arg,
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices_arg) :
+ DynamicModelAC {order_arg, dynamic_g1_sparse_rowval_mx_arg, dynamic_g1_sparse_colval_mx_arg,
+ dynamic_g1_sparse_colptr_mx_arg, dynamic_gN_sparse_indices_arg},
+ DynamicMFilename {modName + ".sparse.dynamic"}
{
}
@@ -61,161 +66,145 @@ DynamicModelMFile::cmplxToReal(mxArray* cmplx_mx)
return real_mx;
}
-void
-DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm)
-{
- int totalCols = mxGetN(sparseMat);
- mwIndex* rowIdxVector = mxGetIr(sparseMat);
- mwIndex* colIdxVector = mxGetJc(sparseMat);
-
- assert(tdm.ncols() == 3);
- /* Under MATLAB, the following check always holds at equality; under Octave,
- there may be an inequality, because Octave diminishes nzmax if one gives
- zeros in the values vector when calling sparse(). */
- assert(tdm.nrows() >= static_cast<int>(mxGetNzmax(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;
- if (!mxIsComplex(sparseMat))
- tdm.get(output_row, 2) = mxGetPr(sparseMat)[rind];
- else
- {
- double real, imag;
-#if MX_HAS_INTERLEAVED_COMPLEX
- mxComplexDouble cmplx = mxGetComplexDoubles(sparseMat)[rind];
- real = cmplx.real;
- imag = cmplx.imag;
-#else
- real = mxGetPr(sparseMat)[rind];
- imag = mxGetPi(sparseMat)[rind];
-#endif
- tdm.get(output_row, 2) = imag == 0.0 ? real : std::numeric_limits<double>::quiet_NaN();
- }
- output_row++;
- }
-
- /* If there are less elements than expected (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 < tdm.nrows())
- {
- 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,
- std::vector<TwoDMatrix>& md) noexcept(false)
+ const std::map<int, int>& dynToDynpp,
+ TensorContainer<FSSparseTensor>& derivatives) noexcept(false)
{
- mxArray* T_m = mxCreateDoubleMatrix(ntt, 1, mxREAL);
+ mxArray* y_mx = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
+ std::copy_n(y.base(), y.length(), mxGetPr(y_mx));
- mxArray* y_m = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
- std::copy_n(y.base(), y.length(), mxGetPr(y_m));
+ mxArray* x_mx = mxCreateDoubleMatrix(1, x.length(), mxREAL);
+ std::copy_n(x.base(), x.length(), mxGetPr(x_mx));
- mxArray* x_m = mxCreateDoubleMatrix(1, x.length(), mxREAL);
- std::copy_n(x.base(), x.length(), mxGetPr(x_m));
+ mxArray* params_mx = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL);
+ std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_mx));
- mxArray* params_m = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL);
- std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_m));
+ mxArray* steady_state_mx = mxCreateDoubleMatrix(ySteady.length(), 1, mxREAL);
+ std::copy_n(ySteady.base(), ySteady.length(), mxGetPr(steady_state_mx));
- mxArray* steady_state_m = mxCreateDoubleMatrix(ySteady.length(), 1, mxREAL);
- std::copy_n(ySteady.base(), ySteady.length(), mxGetPr(steady_state_m));
-
- mxArray* it_m = mxCreateDoubleScalar(1.0);
- mxArray* T_flag_m = mxCreateLogicalScalar(false);
+ mxArray *T_order_mx, *T_mx;
{
- // Compute temporary terms (for all orders)
- std::string funcname = DynamicMFilename + "_g" + std::to_string(md.size()) + "_tt";
- std::array<mxArray*, 1> plhs;
- std::array prhs {T_m, y_m, x_m, params_m, steady_state_m, it_m};
+ // Compute residuals
+ std::string funcname = DynamicMFilename + "_resid";
+ std::array<mxArray*, 3> plhs;
+ std::array prhs {y_mx, x_mx, params_mx, steady_state_mx};
int retVal
= mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
if (retVal != 0)
throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
- mxDestroyArray(T_m);
- T_m = plhs[0];
+ if (!mxIsDouble(plhs[0]) || mxIsSparse(plhs[0]))
+ throw DynareException(__FILE__, __LINE__,
+ "Residual should be a dense array of double floats");
+
+ if (mxIsComplex(plhs[0]))
+ plhs[0] = cmplxToReal(plhs[0]);
+
+ residual = Vector {plhs[0]};
+ mxDestroyArray(plhs[0]);
+
+ T_order_mx = plhs[1];
+ T_mx = plhs[2];
}
{
- // Compute residuals
- std::string funcname = DynamicMFilename + "_resid";
- std::array<mxArray*, 1> plhs;
- std::array prhs {T_m, y_m, x_m, params_m, steady_state_m, it_m, T_flag_m};
+ // Compute Jacobian
+ std::string funcname = DynamicMFilename + "_g1";
+ std::array<mxArray*, 3> plhs;
+ std::array prhs {y_mx,
+ x_mx,
+ params_mx,
+ steady_state_mx,
+ const_cast<mxArray*>(dynamic_g1_sparse_rowval_mx),
+ const_cast<mxArray*>(dynamic_g1_sparse_colval_mx),
+ const_cast<mxArray*>(dynamic_g1_sparse_colptr_mx),
+ T_order_mx,
+ T_mx};
int retVal
= mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
if (retVal != 0)
throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
- if (!mxIsDouble(plhs[0]) || mxIsSparse(plhs[0]))
- throw DynareException(__FILE__, __LINE__,
- "Residual should be a dense array of double floats");
+ assert(static_cast<int>(mxGetN(plhs[0])) == y.length() + x.length());
+ double* g1_v {mxGetPr(plhs[0])};
+ mwIndex* g1_ir {mxGetIr(plhs[0])};
+ mwIndex* g1_jc {mxGetJc(plhs[0])};
- if (mxIsComplex(plhs[0]))
- plhs[0] = cmplxToReal(plhs[0]);
+ IntSequence s(1, 0);
+ auto tensor = std::make_unique<FSSparseTensor>(1, dynToDynpp.size(), residual.length());
+
+ for (size_t j {0}; j < mxGetN(plhs[0]); j++)
+ for (mwIndex k {g1_jc[j]}; k < g1_jc[j + 1]; k++)
+ {
+ s[0] = dynToDynpp.at(j);
+ tensor->insert(s, g1_ir[k], g1_v[k]);
+ }
- residual = Vector {plhs[0]};
mxDestroyArray(plhs[0]);
+
+ mxDestroyArray(T_order_mx);
+ T_order_mx = plhs[1];
+
+ mxDestroyArray(T_mx);
+ T_mx = plhs[2];
+
+ derivatives.insert(std::move(tensor));
}
- for (size_t i = 1; i <= md.size(); i++)
+ for (int o {2}; o <= order; o++)
{
- // Compute model derivatives
- std::string funcname = DynamicMFilename + "_g" + std::to_string(i);
- std::array<mxArray*, 1> plhs;
- std::array prhs {T_m, y_m, x_m, params_m, steady_state_m, it_m, T_flag_m};
+ // Compute higher derivatives
+ std::string funcname = DynamicMFilename + "_g" + std::to_string(o);
+ std::array<mxArray*, 3> plhs;
+ std::array prhs {y_mx, x_mx, params_mx, steady_state_mx, T_order_mx, T_mx};
int retVal
= mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
if (retVal != 0)
throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
- if (!mxIsDouble(plhs[0]))
- throw DynareException(__FILE__, __LINE__,
- "Derivatives matrix at order " + std::to_string(i)
- + "should be an array of double floats");
+ const mxArray* sparse_indices_mx {dynamic_gN_sparse_indices[o - 2]};
+ size_t nnz {mxGetM(sparse_indices_mx)};
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T* sparse_indices {mxGetInt32s(sparse_indices_mx)};
+#else
+ const int32_T* sparse_indices {static_cast<const int32_T*>(mxGetData(sparse_indices_mx))};
+#endif
- if (i == 1)
- {
- if (mxIsSparse(plhs[0]))
- throw DynareException(__FILE__, __LINE__,
- "Derivatives matrix at order " + std::to_string(i)
- + " should be dense");
- assert(static_cast<int>(mxGetM(plhs[0])) == md[i - 1].nrows());
- assert(static_cast<int>(mxGetN(plhs[0])) == md[i - 1].ncols());
- std::copy_n(mxGetPr(plhs[0]), mxGetM(plhs[0]) * mxGetN(plhs[0]), md[i - 1].base());
- }
- else
+ assert(mxGetNumberOfElements(plhs[0]) == nnz);
+ double* gN_v {mxGetPr(plhs[0])};
+
+ IntSequence s(o, 0);
+ auto tensor = std::make_unique<FSSparseTensor>(o, dynToDynpp.size(), residual.length());
+
+ for (size_t k {0}; k < nnz; k++)
{
- if (!mxIsSparse(plhs[0]))
- throw DynareException(__FILE__, __LINE__,
- "Derivatives matrix at order " + std::to_string(i)
- + " should be sparse");
- unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[0], md[i - 1]);
+ for (int i {0}; i < o; i++)
+ s[i] = dynToDynpp.at(sparse_indices[k + (i + 1) * nnz] - 1);
+ s.sort();
+ tensor->insert(s, sparse_indices[k] - 1, gN_v[k]);
}
mxDestroyArray(plhs[0]);
+
+ mxDestroyArray(T_order_mx);
+ T_order_mx = plhs[1];
+
+ mxDestroyArray(T_mx);
+ T_mx = plhs[2];
+
+ derivatives.insert(std::move(tensor));
}
- mxDestroyArray(T_m);
- mxDestroyArray(y_m);
- mxDestroyArray(x_m);
- mxDestroyArray(params_m);
- mxDestroyArray(steady_state_m);
- mxDestroyArray(it_m);
- mxDestroyArray(T_flag_m);
+ mxDestroyArray(y_mx);
+ mxDestroyArray(x_mx);
+ mxDestroyArray(params_mx);
+ mxDestroyArray(steady_state_mx);
+ mxDestroyArray(T_order_mx);
+ mxDestroyArray(T_mx);
}
diff --git a/mex/sources/libkorder/dynamic_m.hh b/mex/sources/libkorder/dynamic_m.hh
index bf89e984758a2906a9f2636a105c2d006331396d..1e00eecdd5dbad6b5a62f3af6e4fd91ff9da12d5 100644
--- a/mex/sources/libkorder/dynamic_m.hh
+++ b/mex/sources/libkorder/dynamic_m.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -22,9 +22,6 @@
#include "dynamic_abstract_class.hh"
-#include "mex.h"
-#include <dynmex.h>
-
/**
* handles calls to <model>_dynamic.m
*
@@ -33,10 +30,6 @@ class DynamicModelMFile : public DynamicModelAC
{
private:
const std::string DynamicMFilename;
- /* Unpack a sparse matrix (of double floats) into a TwoDMatrix object.
- Real elements of the original matrix are copied as-is to the new one.
- Complex elements are replaced by NaNs. */
- static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm);
/* Given a complex dense matrix (of double floats), returns a real dense matrix of same size.
Real elements of the original matrix are copied as-is to the new one.
Complex elements are replaced by NaNs.
@@ -44,8 +37,14 @@ private:
static mxArray* cmplxToReal(mxArray* m);
public:
- explicit DynamicModelMFile(const std::string& modName, int ntt_arg);
+ DynamicModelMFile(const std::string& modName, int order_arg,
+ const mxArray* dynamic_g1_sparse_rowval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colval_mx_arg,
+ const mxArray* dynamic_g1_sparse_colptr_mx_arg,
+ const std::vector<const mxArray*> dynamic_gN_sparse_indices_arg);
void eval(const Vector& y, const Vector& x, const Vector& params, const Vector& ySteady,
- Vector& residual, std::vector<TwoDMatrix>& md) override;
+ Vector& residual, const std::map<int, int>& dynToDynpp,
+ TensorContainer<FSSparseTensor>& derivatives) noexcept(false) override;
};
+
#endif
diff --git a/mex/sources/libkorder/k_ord_dynare.cc b/mex/sources/libkorder/k_ord_dynare.cc
index 2d7d3d2f0d026ff97eb1a29f45b267c01a0e9fab..e1ce47d04cf7f2e1f2bf915f011b951acb70dd18 100644
--- a/mex/sources/libkorder/k_ord_dynare.cc
+++ b/mex/sources/libkorder/k_ord_dynare.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2008-2023 Dynare Team
+ * Copyright © 2008-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -31,7 +31,7 @@ KordpDynare::KordpDynare(const std::vector<std::string>& endo, const std::vector
int nstat, int npred, int nforw, int nboth, const ConstVector& nnzd,
int nsteps, int norder, Journal& jr,
std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
- const std::vector<int>& dr_order, const ConstTwoDMatrix& llincidence) :
+ const std::vector<int>& dr_order) :
nStat {nstat},
nBoth {nboth},
nPred {npred},
@@ -53,7 +53,6 @@ KordpDynare::KordpDynare(const std::vector<std::string>& endo, const std::vector
dnl {endo},
denl {exo},
dsnl {*this, dnl, denl},
- ll_Incidence {llincidence},
dynamicModelFile {std::move(dynamicModelFile_arg)}
{
computeJacobianPermutation(dr_order);
@@ -91,97 +90,17 @@ KordpDynare::calcDerivativesAtSteady()
{
assert(md.begin() == md.end());
- std::vector<TwoDMatrix> dyn_md; // Model derivatives, in Dynare form
-
- dyn_md.emplace_back(nY, nJcols); // Allocate Jacobian
- dyn_md.back().zeros();
-
- for (int i = 2; i <= nOrder; i++)
- {
- // Higher order derivatives, as sparse (3-column) matrices
- dyn_md.emplace_back(static_cast<int>(NNZD[i - 1]), 3);
- dyn_md.back().zeros();
- }
-
Vector xx(nexog());
xx.zeros();
Vector out(nY);
out.zeros();
- Vector llxSteady(nJcols - nExog);
- LLxSteady(ySteady, llxSteady);
-
- dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, dyn_md);
-
- for (int i = 1; i <= nOrder; i++)
- populateDerivativesContainer(dyn_md, i);
-}
-
-void
-KordpDynare::populateDerivativesContainer(const std::vector<TwoDMatrix>& dyn_md, int ord)
-{
- const TwoDMatrix& g = dyn_md[ord - 1];
-
- // model derivatives FSSparseTensor instance
- auto mdTi = std::make_unique<FSSparseTensor>(ord, nJcols, nY);
-
- IntSequence s(ord, 0);
-
- if (ord == 1)
- for (int i = 0; i < g.ncols(); i++)
- {
- for (int j = 0; j < g.nrows(); j++)
- {
- double x = g.get(j, dynppToDyn[s[0]]);
- if (x != 0.0)
- mdTi->insert(s, j, x);
- }
- s[0]++;
- }
- else // ord ≥ 2
- 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;
- if (j < 0 || i1 < 0)
- continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
-
- for (int k = 0; k < ord; k++)
- {
- s[k] = dynToDynpp[i1 % nJcols];
- i1 /= nJcols;
- }
-
- if (ord == 2 && !s.isSorted())
- continue; // Skip symmetric elements (only needed at order 2)
- else if (ord > 2)
- s.sort(); // For higher order, canonicalize the multi-index
-
- double x = g.get(i, 2);
- mdTi->insert(s, j, x);
- }
-
- md.insert(std::move(mdTi));
-}
+ Vector llxSteady(3 * nY);
+ std::copy_n(ySteady.base(), nY, llxSteady.base());
+ std::copy_n(ySteady.base(), nY, llxSteady.base() + nY);
+ std::copy_n(ySteady.base(), nY, llxSteady.base() + 2 * nY);
-/* Returns ySteady extended with leads and lags suitable for passing to
- <model>_dynamic */
-void
-KordpDynare::LLxSteady(const Vector& yS, Vector& llxSteady)
-{
- if (yS.length() == nJcols - nExog)
- throw DynareException(__FILE__, __LINE__, "ySteady already of right size");
-
- /* Create temporary square 2D matrix size nEndo×nEndo (sparse)
- 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))
- llxSteady[static_cast<int>(ll_Incidence.get(ll_row, i)) - 1] = yS[i];
+ dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, dynToDynpp, md);
}
/*
@@ -191,53 +110,33 @@ KordpDynare::LLxSteady(const Vector& yS, Vector& llxSteady)
If one defines:
– y (resp. x) as the vector of all endogenous (size nY), in DR-order (resp.
declaration order)
- – yâ» (resp. xâ») as the vector of endogenous that appear at previous period (size nYs),
- in DR-order (resp. declaration order)
- – y⺠(resp. xâº) as the vector of endogenous that appear at future period (size nYss) in
- DR-order (resp. declaration order)
+ – y⻠as the vector of endogenous that appear at previous period (size nYs),
+ in DR-order
+ – y⺠as the vector of endogenous that appear at future period (size nYss) in
+ DR-order
– u as the vector of exogenous (size nExog)
- In Dynare, the ordering is (xâ», x, xâº, u).
- In Dynare++, the ordering is (yâº, y, yâ», u).
-
- dr_order is typically equal to oo_.dr.order_var.
+ In Dynare++, the vector is of size nY+nYs+nYss+nExog and the ordering is (yâº, y, yâ», u).
+ In Dynare, the vector is of size 3*nY+nExog and the ordering is (x, x, x, u).
*/
void
KordpDynare::computeJacobianPermutation(const std::vector<int>& dr_order)
{
- // Compute restricted inverse DR-orderings: xâ»â†’yâ» and xâºâ†’yâº
- std::vector<int> dr_inv_order_forw(nBoth + nForw), dr_inv_order_pred(nBoth + nPred);
- std::iota(dr_inv_order_forw.begin(), dr_inv_order_forw.end(), 0);
- std::sort(dr_inv_order_forw.begin(), dr_inv_order_forw.end(), [&](int i, int j) {
- return dr_order[nStat + nPred + i] < dr_order[nStat + nPred + j];
- });
- std::iota(dr_inv_order_pred.begin(), dr_inv_order_pred.end(), 0);
- std::sort(dr_inv_order_pred.begin(), dr_inv_order_pred.end(),
- [&](int i, int j) { return dr_order[nStat + i] < dr_order[nStat + j]; });
-
- // Compute restricted DR-orderings: yâ»â†’xâ» and yâºâ†’xâº
- std::vector<int> dr_order_forw(nBoth + nForw), dr_order_pred(nBoth + nPred);
- for (int i = 0; i < nBoth + nForw; i++)
- dr_order_forw[dr_inv_order_forw[i]] = i;
- for (int i = 0; i < nBoth + nPred; i++)
- dr_order_pred[dr_inv_order_pred[i]] = i;
-
// Compute Dynare++ → Dynare ordering
- dynppToDyn.resize(nJcols);
- int j = 0;
+ std::vector<int> dynppToDyn(nJcols);
+ int j {0};
for (; j < nYss; j++)
- dynppToDyn[j] = dr_order_forw[j] + nYs + nY; // Forward variables
+ dynppToDyn[j] = dr_order[j + nStat + nPred] + 2 * nY; // Forward variables
for (; j < nYss + nY; j++)
- dynppToDyn[j] = dr_order[j - nYss] + nYs; // Variables in current period
+ dynppToDyn[j] = dr_order[j - nYss] + nY; // Variables in current period
for (; j < nYss + nY + nYs; j++)
- dynppToDyn[j] = dr_order_pred[j - nY - nYss]; // Predetermined variables
+ dynppToDyn[j] = dr_order[j - (nY + nYss) + nStat]; // Predetermined variables
for (; j < nJcols; j++)
- dynppToDyn[j] = j; // Exogenous
+ dynppToDyn[j] = j - (nY + nYss + nYs) + 3 * nY; // Exogenous
// Compute Dynare → Dynare++ ordering
- dynToDynpp.resize(nJcols);
- for (int i = 0; i < nJcols; i++)
- dynToDynpp[dynppToDyn[i]] = i;
+ for (int i {0}; i < nJcols; i++)
+ dynToDynpp.emplace(dynppToDyn[i], i);
}
DynareNameList::DynareNameList(std::vector<std::string> names_arg) : names(std::move(names_arg))
diff --git a/mex/sources/libkorder/k_ord_dynare.hh b/mex/sources/libkorder/k_ord_dynare.hh
index 6d45ce82c340d1cbd572bc2f6c0a9fd0ca436784..af3324acc0e7f5050e57ef86bcb1d7656a915a90 100644
--- a/mex/sources/libkorder/k_ord_dynare.hh
+++ b/mex/sources/libkorder/k_ord_dynare.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2008-2023 Dynare Team
+ * Copyright © 2008-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -101,9 +101,7 @@ private:
TensorContainer<FSSparseTensor> md; // Model derivatives, in Dynare++ form
DynareNameList dnl, denl;
DynareStateNameList dsnl;
- const ConstTwoDMatrix& ll_Incidence;
- std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
- std::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
+ std::map<int, int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
std::unique_ptr<DynamicModelAC> dynamicModelFile;
@@ -112,7 +110,7 @@ public:
int num_exo, int num_par, Vector& ySteady, TwoDMatrix& vCov, Vector& params,
int nstat, int nPred, int nforw, int nboth, const ConstVector& NNZD, int nSteps,
int ord, Journal& jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
- const std::vector<int>& varOrder, const ConstTwoDMatrix& ll_Incidence);
+ const std::vector<int>& dr_order);
[[nodiscard]] int
nstat() const override
@@ -154,16 +152,6 @@ public:
{
return nOrder;
}
- [[nodiscard]] const std::vector<int>&
- getDynppToDyn() const
- {
- return dynppToDyn;
- }
- [[nodiscard]] const std::vector<int>&
- getDynToDynpp() const
- {
- return dynToDynpp;
- }
[[nodiscard]] const NameList&
getAllEndoNames() const override
{
@@ -214,14 +202,9 @@ public:
}
private:
- // Given the steady state in yS, returns in llxSteady the steady state extended with leads and
- // lags
- 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 model derivatives in Dynare++ form (at a given order) given the Dynare form
- void populateDerivativesContainer(const std::vector<TwoDMatrix>& dyn_md, int ord);
+ void computeJacobianPermutation(const std::vector<int>& dr_order);
};
#endif