diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
index f85f1616e025b08b6c049a125c37c058c1f4b17f..a3148239041b5a2fd3371aed03d5cd625bca7387 100644
--- a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2019-2022 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@@ -22,44 +22,51 @@
#include "DynamicModelCaller.hh"
#include <algorithm>
+#include <filesystem>
+
+using namespace std::literals::string_literals;
std::string DynamicModelCaller::error_msg;
-#if defined(_WIN32) || defined(__CYGWIN32__)
-HINSTANCE DynamicModelDllCaller::dynamic_mex{nullptr};
+#if !defined(_WIN32) && !defined(__CYGWIN32__)
+void *DynamicModelDllCaller::resid_mex{nullptr};
+void *DynamicModelDllCaller::g1_mex{nullptr};
#else
-void *DynamicModelDllCaller::dynamic_mex{nullptr};
+HINSTANCE DynamicModelDllCaller::resid_mex{nullptr};
+HINSTANCE DynamicModelDllCaller::g1_mex{nullptr};
#endif
DynamicModelDllCaller::dynamic_tt_fct DynamicModelDllCaller::residual_tt_fct{nullptr}, DynamicModelDllCaller::g1_tt_fct{nullptr};
-DynamicModelDllCaller::dynamic_deriv_fct DynamicModelDllCaller::residual_fct{nullptr}, DynamicModelDllCaller::g1_fct{nullptr};
+DynamicModelDllCaller::dynamic_fct DynamicModelDllCaller::residual_fct{nullptr}, DynamicModelDllCaller::g1_fct{nullptr};
void
DynamicModelDllCaller::load_dll(const std::string &basename)
{
// Load symbols from dynamic MEX
- std::string mex_name;
-#if !defined(__CYGWIN32__) && !defined(_WIN32)
- mex_name = "./";
-#endif
- mex_name += "+" + basename + "/dynamic" + MEXEXT;
+ const std::filesystem::path sparse_dir {"+" + basename + "/+sparse/"};
+ const std::filesystem::path resid_mex_name {sparse_dir / ("dynamic_resid"s + MEXEXT)},
+ g1_mex_name {sparse_dir / ("dynamic_g1"s + MEXEXT)};
#if !defined(__CYGWIN32__) && !defined(_WIN32)
- dynamic_mex = dlopen(mex_name.c_str(), RTLD_NOW);
+ resid_mex = dlopen(resid_mex_name.c_str(), RTLD_NOW);
+ g1_mex = dlopen(g1_mex_name.c_str(), RTLD_NOW);
#else
- dynamic_mex = LoadLibrary(mex_name.c_str());
+ resid_mex = LoadLibraryW(resid_mex_name.c_str());
+ g1_mex = LoadLibraryW(g1_mex_name.c_str());
#endif
- if (!dynamic_mex)
- mexErrMsgTxt("Can't load dynamic MEX file");
+ if (!resid_mex)
+ mexErrMsgTxt("Can't load dynamic_resid MEX file");
+ if (!g1_mex)
+ mexErrMsgTxt("Can't load dynamic_g1 MEX file");
#if !defined(__CYGWIN32__) && !defined(_WIN32)
- residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_resid_tt"));
- residual_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_resid"));
- g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_g1_tt"));
- g1_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_g1"));
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(resid_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_fct>(dlsym(resid_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(g1_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_fct>(dlsym(g1_mex, "dynamic_g1"));
#else
- residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_resid_tt"));
- residual_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_resid"));
- g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_g1_tt"));
- g1_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_g1"));
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(resid_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_fct>(GetProcAddress(resid_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(g1_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_fct>(GetProcAddress(g1_mex, "dynamic_g1"));
#endif
if (!residual_tt_fct || !residual_fct || !g1_tt_fct || !g1_fct)
mexErrMsgTxt("Can't load functions in dynamic MEX file");
@@ -69,91 +76,120 @@ void
DynamicModelDllCaller::unload_dll()
{
#if !defined(__CYGWIN32__) && !defined(_WIN32)
- dlclose(dynamic_mex);
+ dlclose(resid_mex);
+ dlclose(g1_mex);
#else
- FreeLibrary(dynamic_mex);
+ FreeLibrary(resid_mex);
+ FreeLibrary(g1_mex);
#endif
}
-DynamicModelDllCaller::DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool linear_arg, bool compute_jacobian_arg) :
+DynamicModelDllCaller::DynamicModelDllCaller(size_t ntt, mwIndex ny, mwIndex nx, const double *params_arg, const double *steady_state_arg, const int32_T *g1_sparse_colptr_arg, bool linear_arg, bool compute_jacobian_arg) :
DynamicModelCaller{linear_arg, compute_jacobian_arg},
- nb_row_x{nb_row_x_arg}, x{x_arg}, params{params_arg}, steady_state{steady_state_arg}
+ params{params_arg}, steady_state{steady_state_arg},
+ g1_sparse_colptr{g1_sparse_colptr_arg}
{
tt = std::make_unique<double[]>(ntt);
- y_p = std::make_unique<double[]>(ndynvars);
+ y_p = std::make_unique<double[]>(3*ny);
+ x_p = std::make_unique<double[]>(nx);
if (compute_jacobian)
- jacobian_p = std::make_unique<double[]>((ndynvars+nx)*ny);
+ jacobian_p = std::make_unique<double[]>(g1_sparse_colptr[3*ny+nx]-1);
}
void
-DynamicModelDllCaller::eval(int it, double *resid)
+DynamicModelDllCaller::copy_jacobian_column(mwIndex col, double *dest) const
{
- residual_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
- residual_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), resid);
+ std::copy_n(jacobian_p.get() + g1_sparse_colptr[col]-1,
+ g1_sparse_colptr[col+1] - g1_sparse_colptr[col], dest);
+}
+
+void
+DynamicModelDllCaller::eval(double *resid)
+{
+ residual_tt_fct(y_p.get(), x_p.get(), params, steady_state, tt.get());
+ residual_fct(y_p.get(), x_p.get(), params, steady_state, tt.get(), resid);
if (compute_jacobian)
{
- g1_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
- g1_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), jacobian_p.get());
+ g1_tt_fct(y_p.get(), x_p.get(), params, steady_state, tt.get());
+ g1_fct(y_p.get(), x_p.get(), params, steady_state, tt.get(), jacobian_p.get());
if (linear)
compute_jacobian = false; // If model is linear, no need to recompute Jacobian later
}
}
-DynamicModelMatlabCaller::DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool linear_arg, bool compute_jacobian_arg) :
+DynamicModelMatlabCaller::DynamicModelMatlabCaller(std::string basename_arg, mwIndex ny, mwIndex nx, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, const mxArray *g1_sparse_rowval_mx_arg, const mxArray *g1_sparse_colval_mx_arg, const mxArray *g1_sparse_colptr_mx_arg, bool linear_arg, bool compute_jacobian_arg) :
DynamicModelCaller{linear_arg, compute_jacobian_arg},
basename{std::move(basename_arg)},
- T_mx{mxCreateDoubleMatrix(ntt, 1, mxREAL)},
- y_mx{mxCreateDoubleMatrix(ndynvars, 1, mxREAL)},
- it_mx{mxCreateDoubleScalar(0)},
- T_flag_mx{mxCreateLogicalScalar(false)},
+ y_mx{mxCreateDoubleMatrix(3*ny, 1, mxREAL)},
+ x_mx{mxCreateDoubleMatrix(nx, 1, mxREAL)},
jacobian_mx{nullptr},
- x_mx{mxDuplicateArray(x_mx_arg)},
params_mx{mxDuplicateArray(params_mx_arg)},
- steady_state_mx{mxDuplicateArray(steady_state_mx_arg)}
+ steady_state_mx{mxDuplicateArray(steady_state_mx_arg)},
+ g1_sparse_rowval_mx{mxDuplicateArray(g1_sparse_rowval_mx_arg)},
+ g1_sparse_colval_mx{mxDuplicateArray(g1_sparse_colval_mx_arg)},
+ g1_sparse_colptr_mx{mxDuplicateArray(g1_sparse_colptr_mx_arg)}
{
}
DynamicModelMatlabCaller::~DynamicModelMatlabCaller()
{
- mxDestroyArray(T_mx);
mxDestroyArray(y_mx);
- mxDestroyArray(it_mx);
- mxDestroyArray(T_flag_mx);
+ mxDestroyArray(x_mx);
if (jacobian_mx)
mxDestroyArray(jacobian_mx);
- mxDestroyArray(x_mx);
mxDestroyArray(params_mx);
mxDestroyArray(steady_state_mx);
+ mxDestroyArray(g1_sparse_rowval_mx);
+ mxDestroyArray(g1_sparse_colval_mx);
+ mxDestroyArray(g1_sparse_colptr_mx);
}
void
-DynamicModelMatlabCaller::eval(int it, double *resid)
+DynamicModelMatlabCaller::copy_jacobian_column(mwIndex col, double *dest) const
{
- *mxGetPr(it_mx) = it + 1;
-
- {
- // Compute temporary terms
- std::string funcname = basename + (compute_jacobian ? ".dynamic_g1_tt" : ".dynamic_resid_tt");
- mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx };
+ if (jacobian_mx)
+ {
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T *g1_sparse_rowval {mxGetInt32s(g1_sparse_rowval_mx)};
+ const int32_T *g1_sparse_colptr {mxGetInt32s(g1_sparse_colptr_mx)};
+#else
+ const int32_T *g1_sparse_rowval {static_cast<const int32_T *>(mxGetData(g1_sparse_rowval_mx))};
+ const int32_T *g1_sparse_colptr {static_cast<const int32_T *>(mxGetData(g1_sparse_colptr_mx))};
+#endif
- mxArray *exception = mexCallMATLABWithTrap(1, plhs, 6, prhs, funcname.c_str());
- if (exception)
- {
- error_msg = "An error occurred when calling " + funcname;
- return; // Avoid manipulating null pointers in plhs, see #1832
- }
+ /* We cannot assume that jacobian_mx internally uses
+ g1_sparse_{rowval,colval,colptr}, because the call to sparse() in
+ dynamic_g1.m may have further compressed the matrix by removing
+ elements that are numerically zero, despite being symbolically
+ non-zero. */
+ mwIndex *ir {mxGetIr(jacobian_mx)}, *jc {mxGetJc(jacobian_mx)};
+ mwIndex isrc {jc[col]}; // Index in value array of source Jacobian
+ for (mwIndex idest {0}; // Index in value array of destination Jacobian
+ idest < static_cast<mwIndex>(g1_sparse_colptr[col+1]-g1_sparse_colptr[col]); idest++)
+ {
+ mwIndex row {static_cast<mwIndex>(g1_sparse_rowval[idest+g1_sparse_colptr[col]-1]-1)};
+ while (isrc < jc[col+1] && ir[isrc] < row)
+ isrc++;
+ if (isrc < jc[col+1] && ir[isrc] == row)
+ dest[idest] = mxGetPr(jacobian_mx)[isrc];
+ else
+ dest[idest] = 0.0;
+ }
+ }
+}
- mxDestroyArray(T_mx);
- T_mx = plhs[0];
- }
+void
+DynamicModelMatlabCaller::eval(double *resid)
+{
+ mxArray *T_order_mx, *T_mx;
{
// Compute residuals
- std::string funcname = basename + ".dynamic_resid";
- mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+ std::string funcname {basename + ".sparse.dynamic_resid"};
+ mxArray *plhs[3], *prhs[] = { y_mx, x_mx, params_mx, steady_state_mx };
- mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ mxArray *exception { mexCallMATLABWithTrap(3, plhs, 4, prhs, funcname.c_str()) };
if (exception)
{
error_msg = "An error occurred when calling " + funcname;
@@ -167,19 +203,22 @@ DynamicModelMatlabCaller::eval(int it, double *resid)
}
if (mxIsComplex(plhs[0]))
- plhs[0] = cmplxToReal(plhs[0]);
+ plhs[0] = cmplxToReal<false>(plhs[0]);
std::copy_n(mxGetPr(plhs[0]), mxGetNumberOfElements(plhs[0]), resid);
mxDestroyArray(plhs[0]);
+
+ T_order_mx = plhs[1];
+ T_mx = plhs[2];
}
if (compute_jacobian)
{
// Compute Jacobian
- std::string funcname = basename + ".dynamic_g1";
- mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+ std::string funcname {basename + ".sparse.dynamic_g1"};
+ mxArray *plhs[1], *prhs[] = { y_mx, x_mx, params_mx, steady_state_mx, g1_sparse_rowval_mx, g1_sparse_colval_mx, g1_sparse_colptr_mx, T_order_mx, T_mx };
- mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ mxArray *exception { mexCallMATLABWithTrap(1, plhs, 9, prhs, funcname.c_str()) };
if (exception)
{
error_msg = "An error occurred when calling " + funcname;
@@ -192,48 +231,21 @@ DynamicModelMatlabCaller::eval(int it, double *resid)
jacobian_mx = nullptr;
}
- if (!mxIsDouble(plhs[0]) || mxIsSparse(plhs[0]))
+ if (!mxIsDouble(plhs[0]) || !mxIsSparse(plhs[0]))
{
- error_msg = "Jacobian should be a dense array of double floats";
+ error_msg = "Jacobian should be a sparse array of double floats";
return;
}
if (mxIsComplex(plhs[0]))
- jacobian_mx = cmplxToReal(plhs[0]);
+ jacobian_mx = cmplxToReal<true>(plhs[0]);
else
jacobian_mx = plhs[0];
if (linear)
compute_jacobian = false; // If model is linear, no need to recompute Jacobian later
}
-}
-/* NB: This is a duplicate of DynamicModelMFile::cmplxToReal() in
- k_order_perturbation MEX */
-mxArray *
-DynamicModelMatlabCaller::cmplxToReal(mxArray *cmplx_mx)
-{
- mxArray *real_mx = mxCreateDoubleMatrix(mxGetM(cmplx_mx), mxGetN(cmplx_mx), mxREAL);
-
-#if MX_HAS_INTERLEAVED_COMPLEX
- mxComplexDouble *cmplx = mxGetComplexDoubles(cmplx_mx);
-#else
- double *cmplx_real = mxGetPr(cmplx_mx);
- double *cmplx_imag = mxGetPi(cmplx_mx);
-#endif
- double *real = mxGetPr(real_mx);
-
- for (size_t i = 0; i < mxGetNumberOfElements(cmplx_mx); i++)
-#if MX_HAS_INTERLEAVED_COMPLEX
- if (cmplx[i].imag == 0.0)
- real[i] = cmplx[i].real;
-#else
- if (cmplx_imag[i] == 0.0)
- real[i] = cmplx_real[i];
-#endif
- else
- real[i] = std::numeric_limits<double>::quiet_NaN();
-
- mxDestroyArray(cmplx_mx);
- return real_mx;
+ mxDestroyArray(T_order_mx);
+ mxDestroyArray(T_mx);
}
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
index 2cf682895d3fee4a0deb261d969c21ea6c4a9591..9bbc7246e82c31b4951c53a1840cc9457cb3d578 100644
--- a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2019-2022 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@@ -20,6 +20,8 @@
#include <string>
#include <memory>
#include <limits>
+#include <type_traits>
+#include <algorithm>
#if defined(_WIN32) || defined(__CYGWIN32__)
# ifndef NOMINMAX
@@ -41,43 +43,47 @@ public:
DynamicModelCaller(bool linear_arg, bool compute_jacobian_arg) : linear{linear_arg}, compute_jacobian{compute_jacobian_arg}
{
- };
+ }
virtual ~DynamicModelCaller() = default;
- virtual double &y(size_t i) const = 0;
- virtual double jacobian(size_t i) const = 0;
- virtual void eval(int it, double *resid) = 0;
+ virtual double *y() const = 0;
+ virtual double *x() const = 0;
+ /* Copy a column of the Jacobian to dest.
+ Only copies non-zero elements, according to g1_sparse_{rowval,colval,colptr}. */
+ virtual void copy_jacobian_column(mwIndex col, double *dest) const = 0;
+ virtual void eval(double *resid) = 0;
};
class DynamicModelDllCaller : public DynamicModelCaller
{
private:
#if defined(_WIN32) || defined(__CYGWIN32__)
- static HINSTANCE dynamic_mex;
+ static HINSTANCE resid_mex, g1_mex;
#else
- static void *dynamic_mex;
+ static void *resid_mex, *g1_mex;
#endif
- 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_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 *deriv);
+ 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 *value);
static dynamic_tt_fct residual_tt_fct, g1_tt_fct;
- static dynamic_deriv_fct residual_fct, g1_fct;
- size_t nb_row_x;
- const double *x, *params, *steady_state;
- std::unique_ptr<double[]> tt, y_p, jacobian_p;
+ static dynamic_fct residual_fct, g1_fct;
+ const double *params, *steady_state;
+ std::unique_ptr<double[]> tt, y_p, x_p, jacobian_p;
+ const int32_T *g1_sparse_colptr;
public:
- DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool linear_arg, bool compute_jacobian_arg);
+ DynamicModelDllCaller(size_t ntt, mwIndex ny, mwIndex nx, const double *params_arg, const double *steady_state_arg, const int32_T *g1_sparse_colptr_arg, bool linear_arg, bool compute_jacobian_arg);
virtual ~DynamicModelDllCaller() = default;
- double &
- y(size_t i) const override
+ double *
+ y() const override
{
- return y_p[i];
- };
- double
- jacobian(size_t i) const override
+ return y_p.get();
+ }
+ double *
+ x() const override
{
- return jacobian_p[i];
- };
- void eval(int it, double *resid) override;
+ return x_p.get();
+ }
+ void copy_jacobian_column(mwIndex col, double *dest) const override;
+ void eval(double *resid) override;
static void load_dll(const std::string &basename);
static void unload_dll();
};
@@ -86,32 +92,74 @@ class DynamicModelMatlabCaller : public DynamicModelCaller
{
private:
std::string basename;
- mxArray *T_mx, *y_mx, *it_mx, *T_flag_mx, *jacobian_mx, *x_mx, *params_mx, *steady_state_mx;
- /* Given a complex dense matrix (of double floats), returns a real dense matrix of same size.
+ mxArray *y_mx, *x_mx, *jacobian_mx, *params_mx, *steady_state_mx, *g1_sparse_rowval_mx, *g1_sparse_colval_mx, *g1_sparse_colptr_mx;
+ /* Given a complex matrix (of double floats), returns a sparse matrix of same size.
Real elements of the original matrix are copied as-is to the new one.
Complex elements are replaced by NaNs.
- Destroys the original matrix. */
- static mxArray *cmplxToReal(mxArray *m);
+ Destroys the original matrix.
+ There are two versions, one for dense matrices, another for sparse
+ matrices. */
+ template<bool sparse>
+ static mxArray *cmplxToReal(mxArray *cmplx_mx);
public:
- DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool linear_arg, bool compute_jacobian_arg);
+ DynamicModelMatlabCaller(std::string basename_arg, mwIndex ny, mwIndex nx, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, const mxArray *g1_sparse_rowval_mx_arg, const mxArray *g1_sparse_colval_mx_arg, const mxArray *g1_sparse_colptr_mx_arg, bool linear_arg, bool compute_jacobian_arg);
~DynamicModelMatlabCaller() override;
- double &
- y(size_t i) const override
+ double *
+ y() const override
{
- return mxGetPr(y_mx)[i];
- };
- double
- jacobian(size_t i) const override
+ return mxGetPr(y_mx);
+ }
+ double *
+ x() const override
{
- return jacobian_mx ? mxGetPr(jacobian_mx)[i] : std::numeric_limits<double>::quiet_NaN();
- };
- void eval(int it, double *resid) override;
+ return mxGetPr(x_mx);
+ }
+ void copy_jacobian_column(mwIndex col, double *dest) const override;
+ void eval(double *resid) override;
class Exception
{
public:
const std::string msg;
Exception(std::string msg_arg) : msg{std::move(msg_arg)}
{
- };
+ }
};
};
+
+template<bool sparse>
+mxArray *
+DynamicModelMatlabCaller::cmplxToReal(mxArray *cmplx_mx)
+{
+ mxArray *real_mx {sparse ?
+ mxCreateSparse(mxGetM(cmplx_mx), mxGetN(cmplx_mx), mxGetNzmax(cmplx_mx), mxREAL) :
+ mxCreateDoubleMatrix(mxGetM(cmplx_mx), mxGetN(cmplx_mx), mxREAL)};
+
+ if constexpr(sparse)
+ {
+ std::copy_n(mxGetIr(cmplx_mx), mxGetNzmax(cmplx_mx), mxGetIr(real_mx));
+ std::copy_n(mxGetJc(cmplx_mx), mxGetN(cmplx_mx)+1, mxGetJc(real_mx));
+ }
+
+#if MX_HAS_INTERLEAVED_COMPLEX
+ mxComplexDouble *cmplx {mxGetComplexDoubles(cmplx_mx)};
+#else
+ double *cmplx_real {mxGetPr(cmplx_mx)};
+ double *cmplx_imag {mxGetPi(cmplx_mx)};
+#endif
+ double *real {mxGetPr(real_mx)};
+ for (std::conditional_t<sparse, mwSize, size_t> i {0};
+ i < [&]{ if constexpr(sparse) return mxGetNzmax(cmplx_mx); else return mxGetNumberOfElements(cmplx_mx); }(); // Use a lambda instead of the ternary operator to have the right type (there is no constexpr ternary operator)
+ i++)
+#if MX_HAS_INTERLEAVED_COMPLEX
+ if (cmplx[i].imag == 0.0)
+ real[i] = cmplx[i].real;
+#else
+ if (cmplx_imag[i] == 0.0)
+ real[i] = cmplx_real[i];
+#endif
+ else
+ real[i] = std::numeric_limits<double>::quiet_NaN();
+
+ mxDestroyArray(cmplx_mx);
+ return real_mx;
+}
diff --git a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
index a9bb4bcf299c0d52da56dccca3a7a396c1e2f2f6..ff7eb159c94802be592523e96875c69fc0b47b92 100644
--- a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
+++ b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2019-2022 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@@ -59,23 +59,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt("M_.maximum_lag should be a numeric scalar");
mwIndex maximum_lag = static_cast<mwIndex>(mxGetScalar(maximum_lag_mx));
- const mxArray *maximum_endo_lag_mx = mxGetField(M_mx, 0, "maximum_endo_lag");
- if (!(maximum_endo_lag_mx && mxIsScalar(maximum_endo_lag_mx) && mxIsNumeric(maximum_endo_lag_mx)))
- mexErrMsgTxt("M_.maximum_endo_lag should be a numeric scalar");
- mwIndex maximum_endo_lag = static_cast<mwIndex>(mxGetScalar(maximum_endo_lag_mx));
-
const mxArray *dynamic_tmp_nbr_mx = mxGetField(M_mx, 0, "dynamic_tmp_nbr");
if (!(dynamic_tmp_nbr_mx && mxIsDouble(dynamic_tmp_nbr_mx) && mxGetNumberOfElements(dynamic_tmp_nbr_mx) >= 2)
|| mxIsComplex(dynamic_tmp_nbr_mx) || mxIsSparse(dynamic_tmp_nbr_mx))
mexErrMsgTxt("M_.dynamic_tmp_nbr should be a real dense array of at least 2 elements");
- size_t ntt = mxGetPr(dynamic_tmp_nbr_mx)[0] + mxGetPr(dynamic_tmp_nbr_mx)[1];
-
- 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) == static_cast<size_t>(2+maximum_endo_lag)
- && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(ny))
- || mxIsComplex(lead_lag_incidence_mx) || mxIsSparse(lead_lag_incidence_mx))
- mexErrMsgTxt("M_.lead_lag_incidence should be a real dense matrix with 2+M_.maximum_endo_lag rows and M_.endo_nbr columns");
- const double *lead_lag_incidence = mxGetPr(lead_lag_incidence_mx);
+ size_t ntt {static_cast<size_t>(mxGetPr(dynamic_tmp_nbr_mx)[0]) +
+ (compute_jacobian ? static_cast<size_t>(mxGetPr(dynamic_tmp_nbr_mx)[1]) : 0)};
const mxArray *has_external_function_mx = mxGetField(M_mx, 0, "has_external_function");
if (!(has_external_function_mx && mxIsLogicalScalar(has_external_function_mx)))
@@ -101,41 +90,6 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt("options_.threads.perfect_foresight_problem should be a numeric scalar");
int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
- // Call <model>.dynamic_g1_nz
- mxArray *g1_nz_plhs[3];
- if (mexCallMATLAB(3, g1_nz_plhs, 0, nullptr, (basename + ".dynamic_g1_nz").c_str()) != 0)
- mexErrMsgTxt(("Could not call " + basename + ".dynamic_g1_nz").c_str());
- const mxArray *nzij_pred_mx = g1_nz_plhs[0];
- const mxArray *nzij_current_mx = g1_nz_plhs[1];
- const mxArray *nzij_fwrd_mx = g1_nz_plhs[2];
-
- if (!(mxIsInt32(nzij_pred_mx) && mxGetN(nzij_pred_mx) == 2))
- mexErrMsgTxt("nzij_pred should be an int32 matrix with 2 columns");
- size_t nnz_pred = mxGetM(nzij_pred_mx);
-#if MX_HAS_INTERLEAVED_COMPLEX
- const int32_T *nzij_pred = mxGetInt32s(nzij_pred_mx);
-#else
- const int32_T *nzij_pred = static_cast<const int32_T *>(mxGetData(nzij_pred_mx));
-#endif
-
- if (!(mxIsInt32(nzij_current_mx) && mxGetN(nzij_current_mx) == 2))
- mexErrMsgTxt("nzij_current should be an int32 matrix with 2 columns");
- size_t nnz_current = mxGetM(nzij_current_mx);
-#if MX_HAS_INTERLEAVED_COMPLEX
- const int32_T *nzij_current = mxGetInt32s(nzij_current_mx);
-#else
- const int32_T *nzij_current = static_cast<const int32_T *>(mxGetData(nzij_current_mx));
-#endif
-
- if (!(mxIsInt32(nzij_fwrd_mx) && mxGetN(nzij_fwrd_mx) == 2))
- mexErrMsgTxt("nzij_fwrd should be an int32 matrix with 2 columns");
- size_t nnz_fwrd = mxGetM(nzij_fwrd_mx);
-#if MX_HAS_INTERLEAVED_COMPLEX
- const int32_T *nzij_fwrd = mxGetInt32s(nzij_fwrd_mx);
-#else
- const int32_T *nzij_fwrd = static_cast<const int32_T *>(mxGetData(nzij_fwrd_mx));
-#endif
-
// Check other input and map it to local variables
if (!(mxIsScalar(periods_mx) && mxIsNumeric(periods_mx)))
mexErrMsgTxt("periods should be a numeric scalar");
@@ -159,6 +113,32 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
size_t nb_row_x = mxGetM(exo_path_mx);
const double *exo_path = mxGetPr(exo_path_mx);
+ const mxArray *g1_sparse_rowval_mx {mxGetField(M_mx, 0, "dynamic_g1_sparse_rowval")};
+ if (!(mxIsInt32(g1_sparse_rowval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_rowval should be an int32 vector");
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T *g1_sparse_rowval {mxGetInt32s(g1_sparse_rowval_mx)};
+#else
+ const int32_T *g1_sparse_rowval {static_cast<const int32_T *>(mxGetData(g1_sparse_rowval_mx))};
+#endif
+
+ const mxArray *g1_sparse_colval_mx {mxGetField(M_mx, 0, "dynamic_g1_sparse_colval")};
+ if (!(mxIsInt32(g1_sparse_colval_mx)))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colval should be an int32 vector");
+ if (mxGetNumberOfElements(g1_sparse_colval_mx) != mxGetNumberOfElements(g1_sparse_rowval_mx))
+ mexErrMsgTxt("M_.dynamic_g1_sparse_colval should have the same length as M_.dynamic_g1_sparse_rowval");
+
+ const mxArray *g1_sparse_colptr_mx {mxGetField(M_mx, 0, "dynamic_g1_sparse_colptr")};
+ if (!(mxIsInt32(g1_sparse_colptr_mx) && mxGetNumberOfElements(g1_sparse_colptr_mx) != 3*static_cast<size_t>(ny)+1))
+ mexErrMsgTxt(("M_.dynamic_g1_sparse_colptr should be an int32 vector with " + std::to_string(3*ny+1) + " elements").c_str());
+#if MX_HAS_INTERLEAVED_COMPLEX
+ const int32_T *g1_sparse_colptr {mxGetInt32s(g1_sparse_colptr_mx)};
+#else
+ const int32_T *g1_sparse_colptr {static_cast<const int32_T *>(mxGetData(g1_sparse_colptr_mx))};
+#endif
+ if (static_cast<size_t>(g1_sparse_colptr[3*ny+nx])-1 != mxGetNumberOfElements(g1_sparse_rowval_mx))
+ mexErrMsgTxt("The size of M_.dynamic_g1_sparse_rowval is not consistent with the last element of M_.dynamic_g1_sparse_colptr");
+
if (!(mxIsDouble(params_mx) && mxGetN(params_mx) == 1))
mexErrMsgTxt("params should be a double precision column-vector");
const double *params = mxGetPr(params_mx);
@@ -171,7 +151,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
plhs[0] = mxCreateDoubleMatrix(periods*ny, 1, mxREAL);
double *stacked_residual = mxGetPr(plhs[0]);
- mwIndex nzmax = periods*nnz_current+(periods-1)*(nnz_pred+nnz_fwrd);
+ /* Number of non-zero values in the stacked Jacobian.
+ Contemporaneous derivatives appear at all periods, while lag or lead
+ derivatives appear at all periods except one. */
+ mwSize nzmax {static_cast<mwSize>((g1_sparse_colptr[3*ny]-1)*(periods-1) + (g1_sparse_colptr[2*ny]-g1_sparse_colptr[ny]))};
double *stacked_jacobian = nullptr;
mwIndex *ir = nullptr, *jc = nullptr;
@@ -188,25 +171,20 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mwIndex k = 0;
jc[0] = 0;
for (mwIndex T = 0; T < periods; T++)
- {
- size_t row_pred = 0, row_current = 0, row_fwrd = 0;
- for (int32_T j = 0; j < static_cast<int32_T>(ny); j++)
- {
- if (T != 0)
- while (row_fwrd < nnz_fwrd && nzij_fwrd[row_fwrd+nnz_fwrd]-1 == j)
- ir[k++] = (T-1)*ny + nzij_fwrd[row_fwrd++]-1;
- while (row_current < nnz_current && nzij_current[row_current+nnz_current]-1 == j)
- ir[k++] = T*ny + nzij_current[row_current++]-1;
- if (T != periods-1)
- while (row_pred < nnz_pred && nzij_pred[row_pred+nnz_pred]-1 == j)
- ir[k++] = (T+1)*ny + nzij_pred[row_pred++]-1;
- jc[T*ny+j+1] = k;
- }
- }
+ for (mwIndex j {0}; j < ny; j++) // Column within the period (i.e. variable)
+ {
+ if (T > 0)
+ for (int32_T idx {g1_sparse_colptr[2*ny+j]-1}; idx < g1_sparse_colptr[2*ny+j+1]-1; idx++)
+ ir[k++] = (T-1)*ny + g1_sparse_rowval[idx]-1; // Derivatives w.r.t. y_{t+1} in T-1
+ for (int32_T idx {g1_sparse_colptr[ny+j]-1}; idx < g1_sparse_colptr[ny+j+1]-1; idx++)
+ ir[k++] = T*ny + g1_sparse_rowval[idx]-1; // Derivatives w.r.t. y_t in T
+ if (T < periods-1)
+ for (int32_T idx {g1_sparse_colptr[j]-1}; idx < g1_sparse_colptr[j+1]-1; idx++)
+ ir[k++] = (T+1)*ny + g1_sparse_rowval[idx]-1; // Derivatives w.r.t. y_{t-1} in T+1
+ jc[T*ny+j+1] = k;
+ }
}
- size_t ndynvars = static_cast<size_t>(*std::max_element(lead_lag_incidence, lead_lag_incidence+(maximum_endo_lag+2)*ny));
-
if (use_dll)
DynamicModelDllCaller::load_dll(basename);
@@ -219,55 +197,61 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
// Allocate (thread-private) model evaluator (which allocates space for temporaries)
std::unique_ptr<DynamicModelCaller> m;
if (use_dll)
- m = std::make_unique<DynamicModelDllCaller>(ntt, nx, ny, ndynvars, exo_path, nb_row_x, params, steady_state, linear, compute_jacobian);
+ m = std::make_unique<DynamicModelDllCaller>(ntt, ny, nx, params, steady_state, g1_sparse_colptr, linear, compute_jacobian);
else
- m = std::make_unique<DynamicModelMatlabCaller>(basename, ntt, ndynvars, exo_path_mx, params_mx, steady_state_mx, linear, compute_jacobian);
+ m = std::make_unique<DynamicModelMatlabCaller>(basename, ny, nx, params_mx, steady_state_mx, g1_sparse_rowval_mx, g1_sparse_colval_mx, g1_sparse_colptr_mx, linear, compute_jacobian);
// Main computing loop
#pragma omp for
for (mwIndex T = 0; T < periods; T++)
{
// Fill vector of dynamic variables
- for (mwIndex j = 0; j < maximum_endo_lag+2; j++)
- for (mwIndex i = 0; i < ny; i++)
- {
- int idx = static_cast<int>(lead_lag_incidence[j+i*(2+maximum_endo_lag)])-1;
- if (idx != -1)
- {
- if (T+j == maximum_endo_lag-1)
- m->y(idx) = y0[i];
- else if (T+j == maximum_endo_lag+periods)
- m->y(idx) = yT[i];
- else
- m->y(idx) = y[i+(T+j-maximum_endo_lag)*ny];
- }
- }
+ if (T > 0 && T < periods-1)
+ std::copy_n(y+(T-1)*ny, 3*ny, m->y());
+ else if (T > 0) // Last simulation period
+ {
+ std::copy_n(y+(T-1)*ny, 2*ny, m->y());
+ std::copy_n(yT, ny, m->y() + 2*ny);
+ }
+ else if (T < periods-1) // First simulation period
+ {
+ std::copy_n(y0, ny, m->y());
+ std::copy_n(y+T*ny, 2*ny, m->y() + ny);
+ }
+ else // Special case: periods=1 (and so T=0)
+ {
+ std::copy_n(y0, ny, m->y());
+ std::copy_n(y, ny, m->y() + ny);
+ std::copy_n(yT, ny, m->y() + 2*ny);
+ }
+
+ // Fill exogenous
+ for (mwIndex j {0}; j < nx; j++)
+ m->x()[j] = exo_path[T+maximum_lag + nb_row_x*j];
// Compute the residual and Jacobian, and fill the stacked residual
- m->eval(T+maximum_lag, stacked_residual+T*ny);
+ m->eval(stacked_residual+T*ny);
if (compute_jacobian)
{
// Fill the stacked jacobian
- for (mwIndex col = T > maximum_endo_lag ? (T-maximum_endo_lag)*ny : 0; // We can't use std::max() here, because mwIndex is unsigned under MATLAB
+ for (mwIndex col { T > 1 ? (T-1)*ny : 0 }; // We can't use std::max() here, because mwIndex is unsigned under MATLAB
col < std::min(periods*ny, (T+2)*ny); col++)
{
mwIndex k = jc[col];
while (k < jc[col+1])
{
- if (ir[k] < T*ny)
+ if (ir[k] >= (T+1)*ny)
+ break; // Nothing to copy for this column
+ if (ir[k] >= T*ny)
{
- k++;
- continue;
+ /* Within the current column, this is the first line of
+ the stacked Jacobian that contains elements from the
+ (small) Jacobian just computed, so copy the whole
+ column of the latter to the former. */
+ m->copy_jacobian_column(col - (T-1)*ny, stacked_jacobian + k);
+ break;
}
- if (ir[k] >= (T+1)*ny)
- break;
-
- mwIndex eq = ir[k]-T*ny;
- mwIndex lli_row = col/ny-(T-maximum_endo_lag); // 0, 1 or 2
- mwIndex lli_col = col%ny;
- mwIndex dynvar = static_cast<mwIndex>(lead_lag_incidence[lli_row+lli_col*(2+maximum_endo_lag)])-1;
- stacked_jacobian[k] = m->jacobian(eq+dynvar*ny);
k++;
}
}