From 40ef954802b923dd61a37f6da341641e40d9f6c9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 3 May 2019 19:24:55 +0200
Subject: [PATCH] Add new dynamic_g1_nz.m generated file, for
 perfect_foresight_stacked_solver DLL

---
 src/DynamicModel.cc | 41 +++++++++++++++++++++++++++++++++++++++++
 src/DynamicModel.hh |  4 ++++
 src/ModFile.cc      |  2 ++
 3 files changed, 47 insertions(+)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 421b46c9..81917191 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -2924,6 +2924,47 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
     }
 }
 
+void
+DynamicModel::writeDynamicJacobianNonZeroElts(const string &basename) const
+{
+  vector<pair<int, int>> nzij_pred, nzij_current, nzij_fwrd; // pairs (tsid, equation)
+  for (const auto &it : derivatives[1])
+    {
+      if (symbol_table.getType(getSymbIDByDerivID(it.first[1])) != SymbolType::endogenous)
+        continue;
+      int tsid = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it.first[1]));
+      int lag = getLagByDerivID(it.first[1]);
+      if (lag == -1)
+        nzij_pred.emplace_back(tsid, it.first[0]);
+      else if (lag == 0)
+        nzij_current.emplace_back(tsid, it.first[0]);
+      else
+        nzij_fwrd.emplace_back(tsid, it.first[0]);
+    }
+  sort(nzij_pred.begin(), nzij_pred.end());
+  sort(nzij_current.begin(), nzij_current.end());
+  sort(nzij_fwrd.begin(), nzij_fwrd.end());
+
+  ofstream output{"+" + basename + "/dynamic_g1_nz.m", ios::out | ios::binary};
+  output << "function [nzij_pred, nzij_current, nzij_fwrd] = dynamic_g1_nz()" << endl
+         << "% Returns the coordinates of non-zero elements in the Jacobian, in column-major order, for each lead/lag (only for endogenous)" << endl;
+  auto print_nzij = [&output](const vector<pair<int, int>> &nzij, const string &name)    {
+      output << "  " << name << " = zeros(" << nzij.size() << ", 2, 'int32');" << endl;
+      int idx = 1;
+      for (const auto &it : nzij)
+        {
+          output << "  " << name << "(" << idx << ",1)=" << it.second+1 << ';'
+                 << " " << name << "(" << idx << ",2)=" << it.first+1 << ';' << endl;
+          idx++;
+        }
+    };
+  print_nzij(nzij_pred, "nzij_pred");
+  print_nzij(nzij_current, "nzij_current");
+  print_nzij(nzij_fwrd, "nzij_fwrd");
+  output << "end" << endl;
+  output.close();
+}
+
 void
 DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 83e34e62..324fa8b7 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -365,6 +365,10 @@ public:
   //! Writes file containing parameters derivatives
   void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
+  //! Writes file containing coordinates of non-zero elements in the Jacobian
+  /*! Used by the perfect_foresight_stacked_jacobian MEX */
+  void writeDynamicJacobianNonZeroElts(const string &basename) const;
+
   //! Converts to nonlinear model (only the equations)
   /*! It assumes that the nonlinear model given in argument has just been allocated */
   void toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index b4f85ea0..159bf6b6 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1101,6 +1101,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
           dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mexext, matlabroot, dynareroot, mod_file_struct.order_option, false);
 
           dynamic_model.writeParamsDerivativesFile(basename, false);
+
+          dynamic_model.writeDynamicJacobianNonZeroElts(basename);
         }
 
       // Create steady state file