diff --git a/src/Makefile.am b/src/Makefile.am
index e880905a420a524eaea0c9475b06d975b26e5321..5583a2ae69a10e6a5eddc581f03ff7dc89c1fec8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,8 +42,8 @@ dynare_m_SOURCES = \
 	Statement.hh \
 	ExprNode.cc \
 	ExprNode.hh \
-	MinimumFeedbackSet.cc \
-	MinimumFeedbackSet.hh \
+	VariableDependencyGraph.cc \
+	VariableDependencyGraph.hh \
 	DynareMain.cc \
 	MacroExpandModFile.cc \
 	CodeInterpreter.hh \
diff --git a/src/MinimumFeedbackSet.cc b/src/MinimumFeedbackSet.cc
deleted file mode 100644
index a162ce33b4de92485e9a3dd6c4a8a0ea747fc94f..0000000000000000000000000000000000000000
--- a/src/MinimumFeedbackSet.cc
+++ /dev/null
@@ -1,455 +0,0 @@
-/*
- * Copyright © 2009-2020 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 <iostream>
-
-#include "MinimumFeedbackSet.hh"
-
-using namespace boost;
-
-namespace MFS
-{
-  using color_t = map<boost::graph_traits<AdjacencyList_t>::vertex_descriptor, boost::default_color_type>;
-  using vector_vertex_descriptor_t = vector<AdjacencyList_t::vertex_descriptor>;
-
-  void
-  Suppress(AdjacencyList_t::vertex_descriptor vertex_to_eliminate, AdjacencyList_t &G)
-  {
-    clear_vertex(vertex_to_eliminate, G);
-    remove_vertex(vertex_to_eliminate, G);
-  }
-
-  void
-  Suppress(int vertex_num, AdjacencyList_t &G)
-  {
-    Suppress(vertex(vertex_num, G), G);
-  }
-
-  void
-  Eliminate(AdjacencyList_t::vertex_descriptor vertex_to_eliminate, AdjacencyList_t &G)
-  {
-    if (in_degree(vertex_to_eliminate, G) > 0 && out_degree(vertex_to_eliminate, G) > 0)
-      {
-        AdjacencyList_t::in_edge_iterator it_in, in_end;
-        AdjacencyList_t::out_edge_iterator it_out, out_end;
-        for (tie(it_in, in_end) = in_edges(vertex_to_eliminate, G); it_in != in_end; ++it_in)
-          for (tie(it_out, out_end) = out_edges(vertex_to_eliminate, G); it_out != out_end; ++it_out)
-            {
-              AdjacencyList_t::edge_descriptor ed;
-              bool exist;
-              tie(ed, exist) = edge(source(*it_in, G), target(*it_out, G), G);
-              if (!exist)
-                add_edge(source(*it_in, G), target(*it_out, G), G);
-            }
-      }
-    Suppress(vertex_to_eliminate, G);
-  }
-
-  bool
-  has_cycle_dfs(AdjacencyList_t &g, AdjacencyList_t::vertex_descriptor u, color_t &color, vector<int> &circuit_stack)
-  {
-    property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, g);
-    color[u] = gray_color;
-    graph_traits<AdjacencyList_t>::out_edge_iterator vi, vi_end;
-    for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi)
-      if (color[target(*vi, g)] == white_color && has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
-        {
-          // cycle detected, return immediately
-          circuit_stack.push_back(v_index[target(*vi, g)]);
-          return true;
-        }
-      else if (color[target(*vi, g)] == gray_color)
-        {
-          // *vi is an ancestor!
-          circuit_stack.push_back(v_index[target(*vi, g)]);
-          return true;
-        }
-    color[u] = black_color;
-    return false;
-  }
-
-  bool
-  has_cycle(vector<int> &circuit_stack, AdjacencyList_t &g)
-  {
-    // Initialize color map to white
-    color_t color;
-    graph_traits<AdjacencyList_t>::vertex_iterator vi, vi_end;
-    for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-      color[*vi] = white_color;
-
-    // Perform depth-first search
-    for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-      if (color[*vi] == white_color && has_cycle_dfs(g, *vi, color, circuit_stack))
-        return true;
-
-    return false;
-  }
-
-  void
-  Print(AdjacencyList_t &G)
-  {
-    AdjacencyList_t::vertex_iterator it, it_end;
-    auto v_index = get(vertex_index, G);
-    cout << "Graph\n"
-         << "-----\n";
-    for (tie(it, it_end) = vertices(G); it != it_end; ++it)
-      {
-        cout << "vertex[" << v_index[*it] + 1 << "] <-";
-        AdjacencyList_t::in_edge_iterator it_in, in_end;
-        for (tie(it_in, in_end) = in_edges(*it, G); it_in != in_end; ++it_in)
-          cout << v_index[source(*it_in, G)] + 1 << " ";
-        cout << "\n       ->";
-        AdjacencyList_t::out_edge_iterator it_out, out_end;
-        for (tie(it_out, out_end) = out_edges(*it, G); it_out != out_end; ++it_out)
-          cout << v_index[target(*it_out, G)] + 1 << " ";
-        cout << "\n";
-      }
-  }
-
-  AdjacencyList_t
-  AM_2_AdjacencyList(bool *AM, int n)
-  {
-    AdjacencyList_t G(n);
-    auto v_index = get(vertex_index, G);
-    auto v_index1 = get(vertex_index1, G);
-    for (int i = 0; i < n; i++)
-      {
-        put(v_index, vertex(i, G), i);
-        put(v_index1, vertex(i, G), i);
-      }
-    for (int i = 0; i < n; i++)
-      for (int j = 0; j < n; j++)
-        if (AM[i*n+j])
-          add_edge(vertex(j, G), vertex(i, G), G);
-    return G;
-  }
-
-  AdjacencyList_t
-  extract_subgraph(AdjacencyList_t &G1, set<int> select_index)
-  {
-    int n = select_index.size();
-    AdjacencyList_t G(n);
-    auto v_index = get(vertex_index, G);
-    auto v_index1 = get(vertex_index1, G);
-    auto v1_index = get(vertex_index, G1);
-    map<int, int> reverse_index;
-    set<int>::iterator it;
-    int i;
-    for (it = select_index.begin(), i = 0; i < n; i++, ++it)
-      {
-        reverse_index[get(v1_index, vertex(*it, G1))] = i;
-        put(v_index, vertex(i, G), get(v1_index, vertex(*it, G1)));
-        put(v_index1, vertex(i, G), i);
-      }
-    for (it = select_index.begin(), i = 0; i < n; i++, ++it)
-      {
-        AdjacencyList_t::out_edge_iterator it_out, out_end;
-        auto vi = vertex(*it, G1);
-        for (tie(it_out, out_end) = out_edges(vi, G1); it_out != out_end; ++it_out)
-          {
-            int ii = v1_index[target(*it_out, G1)];
-            if (select_index.find(ii) != select_index.end())
-              add_edge(vertex(reverse_index[get(v1_index, source(*it_out, G1))], G), vertex(reverse_index[get(v1_index, target(*it_out, G1))], G), G);
-          }
-      }
-    return G;
-  }
-
-  vector_vertex_descriptor_t
-  Collect_Doublet(AdjacencyList_t::vertex_descriptor vertex, AdjacencyList_t &G)
-  {
-    AdjacencyList_t::in_edge_iterator it_in, in_end;
-    AdjacencyList_t::out_edge_iterator it_out, out_end;
-    vector<AdjacencyList_t::vertex_descriptor> Doublet;
-    if (in_degree(vertex, G) > 0 && out_degree(vertex, G) > 0)
-      for (tie(it_in, in_end) = in_edges(vertex, G); it_in != in_end; ++it_in)
-        for (tie(it_out, out_end) = out_edges(vertex, G); it_out != out_end; ++it_out)
-          if (source(*it_in, G) == target(*it_out, G) && source(*it_in, G) != target(*it_in, G)) // not a loop
-            Doublet.push_back(source(*it_in, G));
-    return Doublet;
-  }
-
-  bool
-  Vertex_Belong_to_a_Clique(AdjacencyList_t::vertex_descriptor vertex, AdjacencyList_t &G)
-  {
-    vector<AdjacencyList_t::vertex_descriptor> liste;
-    bool agree = true;
-    AdjacencyList_t::in_edge_iterator it_in, in_end;
-    AdjacencyList_t::out_edge_iterator it_out, out_end;
-    tie(it_in, in_end) = in_edges(vertex, G);
-    tie(it_out, out_end) = out_edges(vertex, G);
-    while (it_in != in_end && it_out != out_end && agree)
-      {
-        agree = (source(*it_in, G) == target(*it_out, G) && source(*it_in, G) != target(*it_in, G)); //not a loop
-        liste.push_back(source(*it_in, G));
-        ++it_in;
-        ++it_out;
-      }
-    if (agree)
-      {
-        if (it_in != in_end || it_out != out_end)
-          agree = false;
-        int i = 1;
-        while (i < static_cast<int>(liste.size()) && agree)
-          {
-            int j = i + 1;
-            while (j < static_cast<int>(liste.size()) && agree)
-              {
-                AdjacencyList_t::edge_descriptor ed;
-                bool exist1, exist2;
-                tie(ed, exist1) = edge(liste[i], liste[j], G);
-                tie(ed, exist2) = edge(liste[j], liste[i], G);
-                agree = (exist1 && exist2);
-                j++;
-              }
-            i++;
-          }
-      }
-    return agree;
-  }
-
-  bool
-  Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_t &G)
-  {
-    bool something_has_been_done = false;
-    bool not_a_loop;
-    int i;
-    AdjacencyList_t::vertex_iterator it, it1, ita, it_end;
-    for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
-      {
-        int in_degree_n = in_degree(*it, G);
-        int out_degree_n = out_degree(*it, G);
-        if (in_degree_n <= 1 || out_degree_n <= 1)
-          {
-            not_a_loop = true;
-            if (in_degree_n >= 1 && out_degree_n >= 1) // Do not eliminate a vertex if it loops on itself!
-              {
-                AdjacencyList_t::in_edge_iterator it_in, in_end;
-                for (tie(it_in, in_end) = in_edges(*it, G); it_in != in_end; ++it_in)
-                  if (source(*it_in, G) == target(*it_in, G))
-                    {
-#ifdef verbose
-                      cout << v_index[source(*it_in, G)] << " == " << v_index[target(*it_in, G)] << "\n";
-#endif
-                      not_a_loop = false;
-                    }
-              }
-            if (not_a_loop)
-              {
-#ifdef verbose
-                auto v_index = get(vertex_index, G);
-                cout << "->eliminate vertex[" << v_index[*it] + 1 << "]\n";
-#endif
-                Eliminate(*it, G);
-#ifdef verbose
-                Print(G);
-#endif
-                something_has_been_done = true;
-                if (i > 0)
-                  it = ita;
-                else
-                  {
-                    tie(it, it_end) = vertices(G);
-                    i--;
-                  }
-              }
-          }
-        ita = it;
-      }
-    return something_has_been_done;
-  }
-
-  bool
-  Elimination_of_Vertex_belonging_to_a_clique_Step(AdjacencyList_t &G)
-  {
-    AdjacencyList_t::vertex_iterator it, it1, ita, it_end;
-    bool something_has_been_done = false;
-    int i;
-    for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
-      {
-        if (Vertex_Belong_to_a_Clique(*it, G))
-          {
-#ifdef verbose
-            auto v_index = get(vertex_index, G);
-            cout << "eliminate vertex[" << v_index[*it] + 1 << "]\n";
-#endif
-            Eliminate(*it, G);
-            something_has_been_done = true;
-            if (i > 0)
-              it = ita;
-            else
-              {
-                tie(it, it_end) = vertices(G);
-                i--;
-              }
-          }
-        ita = it;
-      }
-    return something_has_been_done;
-  }
-
-  bool
-  Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(set<int> &feed_back_vertices, AdjacencyList_t &G)
-  {
-    bool something_has_been_done = false;
-    AdjacencyList_t::vertex_iterator it, it_end, ita;
-    int i = 0;
-    for (tie(it, it_end) = vertices(G); it != it_end; ++it, i++)
-      {
-        AdjacencyList_t::edge_descriptor ed;
-        bool exist;
-        tie(ed, exist) = edge(*it, *it, G);
-        if (exist)
-          {
-#ifdef verbose
-            auto v_index = get(vertex_index, G);
-            cout << "store v[*it] = " << v_index[*it]+1 << "\n";
-#endif
-            auto v_index1 = get(vertex_index1, G);
-            feed_back_vertices.insert(v_index1[*it]);
-            Suppress(*it, G);
-            something_has_been_done = true;
-            if (i > 0)
-              it = ita;
-            else
-              {
-                tie(it, it_end) = vertices(G);
-                i--;
-              }
-          }
-        ita = it;
-      }
-    return something_has_been_done;
-  }
-
-  AdjacencyList_t
-  Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices, const AdjacencyList_t &G1)
-  {
-    bool something_has_been_done = true;
-    int cut_ = 0;
-    feed_back_vertices.clear();
-    AdjacencyList_t G(G1);
-    while (num_vertices(G) > 0)
-      {
-        while (something_has_been_done && num_vertices(G) > 0)
-          {
-            //Rule 1
-            something_has_been_done = Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(G) /*or something_has_been_done*/;
-#ifdef verbose
-            cout << "1 something_has_been_done=" << something_has_been_done << "\n";
-#endif
-
-            //Rule 2
-            something_has_been_done = (Elimination_of_Vertex_belonging_to_a_clique_Step(G) || something_has_been_done);
-#ifdef verbose
-            cout << "2 something_has_been_done=" << something_has_been_done << "\n";
-#endif
-
-            //Rule 3
-            something_has_been_done = (Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(feed_back_vertices, G) || something_has_been_done);
-#ifdef verbose
-            cout << "3 something_has_been_done=" << something_has_been_done << "\n";
-#endif
-          }
-        vector<int> circuit;
-        if (!has_cycle(circuit, G))
-          {
-#ifdef verbose
-            cout << "has_cycle=false\n";
-#endif
-            //sort(feed_back_vertices.begin(), feed_back_vertices.end());
-            return G;
-          }
-        if (num_vertices(G) > 0)
-          {
-            /*if nothing has been done in the five previous rule then cut the vertex with the maximum in_degree+out_degree*/
-            int max_degree = 0, num = 0;
-            AdjacencyList_t::vertex_iterator it, it_end, max_degree_index;
-            for (tie(it, it_end) = vertices(G); it != it_end; ++it, num++)
-              if (static_cast<int>(in_degree(*it, G) + out_degree(*it, G)) > max_degree)
-                {
-                  max_degree = in_degree(*it, G) + out_degree(*it, G);
-                  max_degree_index = it;
-                }
-
-            auto v_index1 = get(vertex_index1, G);
-            feed_back_vertices.insert(v_index1[*max_degree_index]);
-            cut_++;
-#ifdef verbose
-            auto v_index = get(vertex_index, G);
-            cout << "--> cut vertex " << v_index[*max_degree_index] + 1 << "\n";
-#endif
-            Suppress(*max_degree_index, G);
-            something_has_been_done = true;
-          }
-      }
-#ifdef verbose
-    cout << "cut_=" << cut_ << "\n";
-#endif
-    //sort(feed_back_vertices.begin(), feed_back_vertices.end());
-    return G;
-  }
-
-  struct rev
-  {
-    bool
-    operator()(const int a, const int b) const
-    {
-      return a > b;
-    }
-  };
-
-  void
-  Reorder_the_recursive_variables(const AdjacencyList_t &G1, set<int> &feedback_vertices, vector< int> &Reordered_Vertices)
-  {
-    AdjacencyList_t G(G1);
-    auto v_index = get(vertex_index, G);
-    set<int, rev> fv;
-    for (auto its = feedback_vertices.begin(); its != feedback_vertices.end(); ++its)
-      fv.insert(*its);
-    int i = 0;
-    for (auto its = fv.begin(); its != fv.end(); ++its, i++)
-      Suppress(*its, G);
-    bool something_has_been_done = true;
-    while (something_has_been_done)
-      {
-        something_has_been_done = false;
-        AdjacencyList_t::vertex_iterator it, it_end, ita;
-        for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
-          {
-            if (in_degree(*it, G) == 0)
-              {
-                Reordered_Vertices.push_back(v_index[*it]);
-                Suppress(*it, G);
-                something_has_been_done = true;
-                if (i > 0)
-                  it = ita;
-                else
-                  {
-                    tie(it, it_end) = vertices(G);
-                    i--;
-                  }
-              }
-            ita = it;
-          }
-      }
-    if (num_vertices(G))
-      cout << "Error in the computation of feedback vertex set\n";
-  }
-}
diff --git a/src/MinimumFeedbackSet.hh b/src/MinimumFeedbackSet.hh
deleted file mode 100644
index 9e10ab86cf20f5b97a3ae8b37165e07ea9400a0f..0000000000000000000000000000000000000000
--- a/src/MinimumFeedbackSet.hh
+++ /dev/null
@@ -1,60 +0,0 @@
-/*
- * Copyright © 2009-2020 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/>.
- */
-
-#ifndef _MINIMUMFEEDBACKSET_HH
-#define _MINIMUMFEEDBACKSET_HH
-
-#include <map>
-#include <vector>
-
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wold-style-cast"
-#include <boost/graph/adjacency_list.hpp>
-#pragma GCC diagnostic pop
-
-using namespace std;
-
-namespace MFS
-{
-  using VertexProperty_t = boost::property<boost::vertex_index_t, int,
-                                           boost::property<boost::vertex_index1_t, int,
-                                                           boost::property<boost::vertex_degree_t, int,
-                                                                           boost::property<boost::vertex_in_degree_t, int,
-                                                                                           boost::property<boost::vertex_out_degree_t, int >>>>>;
-  using AdjacencyList_t = boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>;
-
-  //! Extracts a subgraph
-  /*!
-    \param[in] G1 The original graph
-    \param[in] select_index The vertex indices to select
-    \return The subgraph
-
-    The property vertex_index of the subgraph contains indices of the original
-    graph, the property vertex_index1 contains new contiguous indices specific
-    to the subgraph.
-  */
-  AdjacencyList_t extract_subgraph(AdjacencyList_t &G1, set<int> select_index);
-  //! Return the feedback set
-  AdjacencyList_t Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices, const AdjacencyList_t &G);
-  //! Reorder the recursive variables
-  /*! They appear first in a quasi triangular form and they are followed by the feedback variables */
-  void Reorder_the_recursive_variables(const AdjacencyList_t &G1, set<int> &feedback_vertices, vector<int> &Reordered_Vertices);
-};
-
-#endif // _MINIMUMFEEDBACKSET_HH
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 769665f950f7c720dd3d5b587803a977477ec96b..fc8f899615f4b09def3547120e59fef2a8b74bb7 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -18,14 +18,13 @@
  */
 
 #include "ModelTree.hh"
-#include "MinimumFeedbackSet.hh"
+#include "VariableDependencyGraph.hh"
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wold-style-cast"
 #pragma GCC diagnostic ignored "-Wsign-compare"
 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/max_cardinality_matching.hpp>
-#include <boost/graph/strong_components.hpp>
 #include <boost/graph/topological_sort.hpp>
 #pragma GCC diagnostic pop
 
@@ -740,93 +739,56 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
   int nb_var = variable_reordered.size();
   int n = nb_var - prologue - epilogue;
 
-  MFS::AdjacencyList_t G2(n);
-
-  // It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container
-  auto v_index = get(boost::vertex_index, G2);
-  for (int i = 0; i < n; i++)
-    put(v_index, vertex(i, G2), i);
+  /* Construct the graph representing the dependencies between all
+     variables that do not belong to the prologue or the epilogue. */
+  VariableDependencyGraph G(n);
 
   vector<int> reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var);
-
   for (int i = 0; i < nb_var; i++)
     {
       reverse_equation_reordered[equation_reordered[i]] = i;
       reverse_variable_reordered[variable_reordered[i]] = i;
     }
+
   jacob_map_t tmp_normalized_contemporaneous_jacobian;
   if (cutoff == 0)
-    {
-      set<pair<int, int>> endo;
-      for (int i = 0; i < nb_var; i++)
-        {
-          endo.clear();
-          equations[i]->collectEndogenous(endo);
-          for (const auto &it : endo)
-            tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
-        }
-    }
+    for (int i = 0; i < nb_var; i++)
+      {
+        set<pair<int, int>> endo;
+        equations[i]->collectEndogenous(endo);
+        for (const auto &it : endo)
+          tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
+      }
   else
     tmp_normalized_contemporaneous_jacobian = static_jacobian;
   for (const auto &[key, value] : tmp_normalized_contemporaneous_jacobian)
     if (reverse_equation_reordered[key.first] >= prologue && reverse_equation_reordered[key.first] < nb_var - epilogue
         && reverse_variable_reordered[key.second] >= prologue && reverse_variable_reordered[key.second] < nb_var - epilogue
         && key.first != endo2eq[key.second])
-      add_edge(vertex(reverse_equation_reordered[endo2eq[key.second]]-prologue, G2),
-               vertex(reverse_equation_reordered[key.first]-prologue, G2),
-               G2);
+      add_edge(vertex(reverse_equation_reordered[endo2eq[key.second]]-prologue, G),
+               vertex(reverse_equation_reordered[key.first]-prologue, G),
+               G);
 
-  vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
-  boost::iterator_property_map<int *, boost::property_map<MFS::AdjacencyList_t, boost::vertex_index_t>::type, int, int &> endo2block_map(&endo2block[0], get(boost::vertex_index, G2));
+  /* Compute the mapping between endogenous and blocks, using a strongly
+     connected components (SCC) decomposition */
+  auto [num_scc, endo2block] = G.sortedStronglyConnectedComponents();
 
-  // Compute strongly connected components
-  int num = strong_components(G2, endo2block_map);
+  vector<pair<int, int>> blocks(num_scc, { 0, 0 });
 
-  vector<pair<int, int>> blocks(num, { 0, 0 });
-
-  // Create directed acyclic graph associated to the strongly connected components
-  using DirectedGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS>;
-  DirectedGraph dag(num);
-
-  for (int i = 0; i < static_cast<int>(num_vertices(G2)); i++)
-    {
-      MFS::AdjacencyList_t::out_edge_iterator it_out, out_end;
-      MFS::AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
-      for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
-        {
-          int t_b = endo2block_map[target(*it_out, G2)];
-          int s_b = endo2block_map[source(*it_out, G2)];
-          if (s_b != t_b)
-            add_edge(s_b, t_b, dag);
-        }
-    }
-
-  // Compute topological sort of DAG (ordered list of unordered SCC)
-  deque<int> ordered2unordered;
-  topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
-
-  // Construct mapping from unordered SCC to ordered SCC
-  vector<int> unordered2ordered(num);
-  for (int i = 0; i < num; i++)
-    unordered2ordered[ordered2unordered[i]] = i;
-
-  //This vector contains for each block:
-  //   - first set = equations belonging to the block,
-  //   - second set = the feeback variables,
-  //   - third vector = the reordered non-feedback variables.
-  vector<tuple<set<int>, set<int>, vector<int>>> components_set(num);
+  /* This vector contains for each block:
+     – first set = equations belonging to the block,
+     — second set = the feeback variables,
+     — third vector = the reordered non-feedback variables. */
+  vector<tuple<set<int>, set<int>, vector<int>>> components_set(num_scc);
   for (int i = 0; i < static_cast<int>(endo2block.size()); i++)
     {
-      endo2block[i] = unordered2ordered[endo2block[i]];
       blocks[endo2block[i]].first++;
       get<0>(components_set[endo2block[i]]).insert(i);
     }
 
-  auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block, num);
+  auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block, num_scc);
 
-  vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
-  int order = prologue;
-  //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
+  // Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
   if (select_feedback_variable)
     {
       for (int i = 0; i < n; i++)
@@ -836,17 +798,18 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
             || equation_lag_lead[equation_reordered[i+prologue]].second > 0
             || equation_lag_lead[equation_reordered[i+prologue]].first > 0
             || mfs == 0)
-          add_edge(vertex(i, G2), vertex(i, G2), G2);
+          add_edge(vertex(i, G), vertex(i, G), G);
     }
   else
     for (int i = 0; i < n; i++)
       if (Equation_Type[equation_reordered[i+prologue]].first == EquationType::solve || mfs == 0)
-        add_edge(vertex(i, G2), vertex(i, G2), G2);
+        add_edge(vertex(i, G), vertex(i, G), G);
 
-  //Determines the dynamic structure of each equation
-  vector<int> n_static(prologue+num+epilogue, 0), n_forward(prologue+num+epilogue, 0),
-    n_backward(prologue+num+epilogue, 0), n_mixed(prologue+num+epilogue, 0);
+  // Determines the dynamic structure of each equation
+  vector<int> n_static(prologue+num_scc+epilogue, 0), n_forward(prologue+num_scc+epilogue, 0),
+    n_backward(prologue+num_scc+epilogue, 0), n_mixed(prologue+num_scc+epilogue, 0);
 
+  vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
   for (int i = 0; i < prologue; i++)
     if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
       n_mixed[i]++;
@@ -857,24 +820,23 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
     else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
       n_static[i]++;
 
-  //For each block, the minimum set of feedback variable is computed
-  // and the non-feedback variables are reordered to get
-  // a sub-recursive block without feedback variables
+  /* For each block, the minimum set of feedback variable is computed and the
+     non-feedback variables are reordered to get a sub-recursive block without
+     feedback variables. */
 
-  for (int i = 0; i < num; i++)
+  int order = prologue;
+  for (int i = 0; i < num_scc; i++)
     {
-      MFS::AdjacencyList_t G = MFS::extract_subgraph(G2, get<0>(components_set[i]));
-      set<int> feed_back_vertices;
-      MFS::AdjacencyList_t G1 = MFS::Minimal_set_of_feedback_vertex(feed_back_vertices, G);
-      auto v_index = get(boost::vertex_index, G);
+      auto subG = G.extractSubgraph(get<0>(components_set[i]));
+      auto [G1, feed_back_vertices] = subG.minimalSetOfFeedbackVertices();
+      auto v_index1 = get(boost::vertex_index1, subG);
       get<1>(components_set[i]) = feed_back_vertices;
       blocks[i].second = feed_back_vertices.size();
-      vector<int> Reordered_Vertice;
-      MFS::Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
+      auto reordered_vertices = subG.reorderRecursiveVariables(feed_back_vertices);
 
-      //First we have the recursive equations conditional on feedback variables
+      // First we have the recursive equations conditional on feedback variables
       for (int j = 0; j < 4; j++)
-        for (int its : Reordered_Vertice)
+        for (int its : reordered_vertices)
           {
             bool something_done = false;
             if (j == 2 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second != 0)
@@ -905,36 +867,36 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
               }
           }
 
-      get<2>(components_set[i]) = Reordered_Vertice;
-      //Second we have the equations related to the feedback variables
+      get<2>(components_set[i]) = reordered_vertices;
+      // Second we have the equations related to the feedback variables
       for (int j = 0; j < 4; j++)
         for (int feed_back_vertice : feed_back_vertices)
           {
             bool something_done = false;
-            if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
+            if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second != 0)
               {
                 n_mixed[prologue+i]++;
                 something_done = true;
               }
-            else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
+            else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second != 0)
               {
                 n_forward[prologue+i]++;
                 something_done = true;
               }
-            else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
+            else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second == 0)
               {
                 n_backward[prologue+i]++;
                 something_done = true;
               }
-            else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
+            else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second == 0)
               {
                 n_static[prologue+i]++;
                 something_done = true;
               }
             if (something_done)
               {
-                equation_reordered[order] = tmp_equation_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
-                variable_reordered[order] = tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
+                equation_reordered[order] = tmp_equation_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue];
+                variable_reordered[order] = tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue];
                 order++;
               }
           }
@@ -942,13 +904,13 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
 
   for (int i = 0; i < epilogue; i++)
     if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
-      n_mixed[prologue+num+i]++;
+      n_mixed[prologue+num_scc+i]++;
     else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
-      n_forward[prologue+num+i]++;
+      n_forward[prologue+num_scc+i]++;
     else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
-      n_backward[prologue+num+i]++;
+      n_backward[prologue+num_scc+i]++;
     else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
-      n_static[prologue+num+i]++;
+      n_static[prologue+num_scc+i]++;
 
   inv_equation_reordered.resize(nb_var);
   inv_variable_reordered.resize(nb_var);
diff --git a/src/VariableDependencyGraph.cc b/src/VariableDependencyGraph.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f0187704d6c349f38df01c1af49ffb163be135c6
--- /dev/null
+++ b/src/VariableDependencyGraph.cc
@@ -0,0 +1,468 @@
+/*
+ * Copyright © 2009-2020 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 <iostream>
+#include <algorithm>
+
+#include "VariableDependencyGraph.hh"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wold-style-cast"
+#pragma GCC diagnostic ignored "-Wsign-compare"
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#include <boost/graph/strong_components.hpp>
+#include <boost/graph/topological_sort.hpp>
+#pragma GCC diagnostic pop
+
+using namespace boost;
+
+VariableDependencyGraph::VariableDependencyGraph(int n) : base(n)
+{
+  /* It is necessary to manually initialize the vertex_index property since
+     this graph uses listS and not vecS as underlying vertex container */
+  auto v_index = get(vertex_index, *this);
+  for (int i = 0; i < n; i++)
+    put(v_index, vertex(i, *this), i);
+}
+
+void
+VariableDependencyGraph::suppress(vertex_descriptor vertex_to_eliminate)
+{
+  clear_vertex(vertex_to_eliminate, *this);
+  remove_vertex(vertex_to_eliminate, *this);
+}
+
+void
+VariableDependencyGraph::suppress(int vertex_num)
+{
+  suppress(vertex(vertex_num, *this));
+}
+
+void
+VariableDependencyGraph::eliminate(vertex_descriptor vertex_to_eliminate)
+{
+  if (in_degree(vertex_to_eliminate, *this) > 0 && out_degree(vertex_to_eliminate, *this) > 0)
+    {
+      in_edge_iterator it_in, in_end;
+      out_edge_iterator it_out, out_end;
+      for (tie(it_in, in_end) = in_edges(vertex_to_eliminate, *this); it_in != in_end; ++it_in)
+        for (tie(it_out, out_end) = out_edges(vertex_to_eliminate, *this); it_out != out_end; ++it_out)
+          {
+            edge_descriptor ed;
+            bool exist;
+            tie(ed, exist) = edge(source(*it_in, *this), target(*it_out, *this), *this);
+            if (!exist)
+              add_edge(source(*it_in, *this), target(*it_out, *this), *this);
+          }
+    }
+  suppress(vertex_to_eliminate);
+}
+
+bool
+VariableDependencyGraph::hasCycleDFS(vertex_descriptor u, color_t &color, vector<int> &circuit_stack) const
+{
+  auto v_index = get(vertex_index, *this);
+  color[u] = gray_color;
+  out_edge_iterator vi, vi_end;
+  for (tie(vi, vi_end) = out_edges(u, *this); vi != vi_end; ++vi)
+    if (color[target(*vi, *this)] == white_color && hasCycleDFS(target(*vi, *this), color, circuit_stack))
+      {
+        // cycle detected, return immediately
+        circuit_stack.push_back(v_index[target(*vi, *this)]);
+        return true;
+      }
+    else if (color[target(*vi, *this)] == gray_color)
+      {
+        // *vi is an ancestor!
+        circuit_stack.push_back(v_index[target(*vi, *this)]);
+        return true;
+      }
+  color[u] = black_color;
+  return false;
+}
+
+bool
+VariableDependencyGraph::hasCycle() const
+{
+  // Initialize color map to white
+  color_t color;
+  vector<int> circuit_stack;
+  vertex_iterator vi, vi_end;
+  for (tie(vi, vi_end) = vertices(*this); vi != vi_end; ++vi)
+    color[*vi] = white_color;
+
+  // Perform depth-first search
+  for (tie(vi, vi_end) = vertices(*this); vi != vi_end; ++vi)
+    if (color[*vi] == white_color && hasCycleDFS(*vi, color, circuit_stack))
+      return true;
+
+  return false;
+}
+
+void
+VariableDependencyGraph::print() const
+{
+  vertex_iterator it, it_end;
+  auto v_index = get(vertex_index, *this);
+  cout << "Graph\n"
+       << "-----\n";
+  for (tie(it, it_end) = vertices(*this); it != it_end; ++it)
+    {
+      cout << "vertex[" << v_index[*it] + 1 << "] <-";
+      in_edge_iterator it_in, in_end;
+      for (tie(it_in, in_end) = in_edges(*it, *this); it_in != in_end; ++it_in)
+        cout << v_index[source(*it_in, *this)] + 1 << " ";
+      cout << "\n       ->";
+      out_edge_iterator it_out, out_end;
+      for (tie(it_out, out_end) = out_edges(*it, *this); it_out != out_end; ++it_out)
+        cout << v_index[target(*it_out, *this)] + 1 << " ";
+      cout << "\n";
+    }
+}
+
+VariableDependencyGraph
+VariableDependencyGraph::extractSubgraph(const set<int> &select_index) const
+{
+  int n = select_index.size();
+  VariableDependencyGraph G(n);
+  auto v_index = get(vertex_index, *this);
+  auto v_index1_G = get(vertex_index1, G); // Maps new vertices to original indices
+  map<int, int> reverse_index; // Maps orig indices to new ones
+  set<int>::iterator it;
+  int i;
+  for (it = select_index.begin(), i = 0; i < n; i++, ++it)
+    {
+      reverse_index[get(v_index, vertex(*it, *this))] = i;
+      put(v_index1_G, vertex(i, G), get(v_index, vertex(*it, *this)));
+    }
+  for (it = select_index.begin(), i = 0; i < n; i++, ++it)
+    {
+      out_edge_iterator it_out, out_end;
+      auto vi = vertex(*it, *this);
+      for (tie(it_out, out_end) = out_edges(vi, *this); it_out != out_end; ++it_out)
+        {
+          int ii = v_index[target(*it_out, *this)];
+          if (select_index.find(ii) != select_index.end())
+            add_edge(vertex(reverse_index[get(v_index, source(*it_out, *this))], G),
+                     vertex(reverse_index[get(v_index, target(*it_out, *this))], G), G);
+        }
+    }
+  return G;
+}
+
+bool
+VariableDependencyGraph::vertexBelongsToAClique(vertex_descriptor vertex) const
+{
+  vector<vertex_descriptor> liste;
+  bool agree = true;
+  in_edge_iterator it_in, in_end;
+  out_edge_iterator it_out, out_end;
+  tie(it_in, in_end) = in_edges(vertex, *this);
+  tie(it_out, out_end) = out_edges(vertex, *this);
+  while (it_in != in_end && it_out != out_end && agree)
+    {
+      agree = (source(*it_in, *this) == target(*it_out, *this)
+               && source(*it_in, *this) != target(*it_in, *this)); //not a loop
+      liste.push_back(source(*it_in, *this));
+      ++it_in;
+      ++it_out;
+    }
+  if (agree)
+    {
+      if (it_in != in_end || it_out != out_end)
+        agree = false;
+      int i = 1;
+      while (i < static_cast<int>(liste.size()) && agree)
+        {
+          int j = i + 1;
+          while (j < static_cast<int>(liste.size()) && agree)
+            {
+              edge_descriptor ed;
+              bool exist1, exist2;
+              tie(ed, exist1) = edge(liste[i], liste[j], *this);
+              tie(ed, exist2) = edge(liste[j], liste[i], *this);
+              agree = (exist1 && exist2);
+              j++;
+            }
+          i++;
+        }
+    }
+  return agree;
+}
+
+bool
+VariableDependencyGraph::eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree()
+{
+  bool something_has_been_done = false;
+  bool not_a_loop;
+  int i;
+  vertex_iterator it, it1, ita, it_end;
+  for (tie(it, it_end) = vertices(*this), i = 0; it != it_end; ++it, i++)
+    {
+      int in_degree_n = in_degree(*it, *this);
+      int out_degree_n = out_degree(*it, *this);
+      if (in_degree_n <= 1 || out_degree_n <= 1)
+        {
+          not_a_loop = true;
+          if (in_degree_n >= 1 && out_degree_n >= 1) // Do not eliminate a vertex if it loops on itself!
+            {
+              in_edge_iterator it_in, in_end;
+              for (tie(it_in, in_end) = in_edges(*it, *this); it_in != in_end; ++it_in)
+                if (source(*it_in, *this) == target(*it_in, *this))
+                  {
+#ifdef verbose
+                    cout << v_index[source(*it_in, *this)] << " == " << v_index[target(*it_in, *this)] << "\n";
+#endif
+                    not_a_loop = false;
+                  }
+            }
+          if (not_a_loop)
+            {
+#ifdef verbose
+              auto v_index1 = get(vertex_index1, *this);
+              cout << "->eliminate vertex[" << v_index1[*it] + 1 << "]\n";
+#endif
+              eliminate(*it);
+#ifdef verbose
+              print();
+#endif
+              something_has_been_done = true;
+              if (i > 0)
+                it = ita;
+              else
+                {
+                  tie(it, it_end) = vertices(*this);
+                  i--;
+                }
+            }
+        }
+      ita = it;
+    }
+  return something_has_been_done;
+}
+
+bool
+VariableDependencyGraph::eliminationOfVerticesBelongingToAClique()
+{
+  vertex_iterator it, it1, ita, it_end;
+  bool something_has_been_done = false;
+  int i;
+  for (tie(it, it_end) = vertices(*this), i = 0; it != it_end; ++it, i++)
+    {
+      if (vertexBelongsToAClique(*it))
+        {
+#ifdef verbose
+          auto v_index1 = get(vertex_index1, *this);
+          cout << "eliminate vertex[" << v_index1[*it] + 1 << "]\n";
+#endif
+          eliminate(*it);
+          something_has_been_done = true;
+          if (i > 0)
+            it = ita;
+          else
+            {
+              tie(it, it_end) = vertices(*this);
+              i--;
+            }
+        }
+      ita = it;
+    }
+  return something_has_been_done;
+}
+
+bool
+VariableDependencyGraph::suppressionOfVerticesWithLoop(set<int> &feed_back_vertices)
+{
+  bool something_has_been_done = false;
+  vertex_iterator it, it_end, ita;
+  int i = 0;
+  for (tie(it, it_end) = vertices(*this); it != it_end; ++it, i++)
+    {
+      edge_descriptor ed;
+      bool exist;
+      tie(ed, exist) = edge(*it, *it, *this);
+      if (exist)
+        {
+#ifdef verbose
+          auto v_index1 = get(vertex_index1, *this);
+          cout << "store v[*it] = " << v_index1[*it]+1 << "\n";
+#endif
+          auto v_index = get(vertex_index, *this);
+          feed_back_vertices.insert(v_index[*it]);
+          suppress(*it);
+          something_has_been_done = true;
+          if (i > 0)
+            it = ita;
+          else
+            {
+              tie(it, it_end) = vertices(*this);
+              i--;
+            }
+        }
+      ita = it;
+    }
+  return something_has_been_done;
+}
+
+pair<VariableDependencyGraph, set<int>>
+VariableDependencyGraph::minimalSetOfFeedbackVertices() const
+{
+  bool something_has_been_done = true;
+  int cut_ = 0;
+  set<int> feed_back_vertices;
+  VariableDependencyGraph G(*this);
+  while (num_vertices(G) > 0)
+    {
+      while (something_has_been_done && num_vertices(G) > 0)
+        {
+          //Rule 1
+          something_has_been_done = G.eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree() /*or something_has_been_done*/;
+#ifdef verbose
+          cout << "1 something_has_been_done=" << something_has_been_done << "\n";
+#endif
+
+          //Rule 2
+          something_has_been_done = (G.eliminationOfVerticesBelongingToAClique() || something_has_been_done);
+#ifdef verbose
+          cout << "2 something_has_been_done=" << something_has_been_done << "\n";
+#endif
+
+          //Rule 3
+          something_has_been_done = (G.suppressionOfVerticesWithLoop(feed_back_vertices) || something_has_been_done);
+#ifdef verbose
+          cout << "3 something_has_been_done=" << something_has_been_done << "\n";
+#endif
+        }
+      if (!G.hasCycle())
+        {
+#ifdef verbose
+          cout << "has_cycle=false\n";
+#endif
+          return { G, feed_back_vertices };
+        }
+      if (num_vertices(G) > 0)
+        {
+          /* If nothing has been done in the five previous rule then cut the
+             vertex with the maximum in_degree+out_degree */
+          int max_degree = 0, num = 0;
+          vertex_iterator it, it_end, max_degree_index;
+          for (tie(it, it_end) = vertices(G); it != it_end; ++it, num++)
+            if (static_cast<int>(in_degree(*it, G) + out_degree(*it, G)) > max_degree)
+              {
+                max_degree = in_degree(*it, G) + out_degree(*it, G);
+                max_degree_index = it;
+              }
+
+          auto v_index = get(vertex_index, G);
+          feed_back_vertices.insert(v_index[*max_degree_index]);
+          cut_++;
+#ifdef verbose
+          auto v_index1 = get(vertex_index1, G);
+          cout << "--> cut vertex " << v_index1[*max_degree_index] + 1 << "\n";
+#endif
+          G.suppress(*max_degree_index);
+          something_has_been_done = true;
+        }
+    }
+#ifdef verbose
+  cout << "cut_=" << cut_ << "\n";
+#endif
+  return { G, feed_back_vertices };
+}
+
+vector<int>
+VariableDependencyGraph::reorderRecursiveVariables(const set<int> &feedback_vertices) const
+{
+  vector<int> reordered_vertices;
+  VariableDependencyGraph G(*this);
+  auto v_index = get(vertex_index, G);
+  set<int, greater<int>> fv;
+  for (auto its = feedback_vertices.begin(); its != feedback_vertices.end(); ++its)
+    fv.insert(*its);
+  int i = 0;
+  for (auto its = fv.begin(); its != fv.end(); ++its, i++)
+    G.suppress(*its);
+  bool something_has_been_done = true;
+  while (something_has_been_done)
+    {
+      something_has_been_done = false;
+      vertex_iterator it, it_end, ita;
+      for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
+        {
+          if (in_degree(*it, G) == 0)
+            {
+              reordered_vertices.push_back(v_index[*it]);
+              G.suppress(*it);
+              something_has_been_done = true;
+              if (i > 0)
+                it = ita;
+              else
+                {
+                  tie(it, it_end) = vertices(G);
+                  i--;
+                }
+            }
+          ita = it;
+        }
+    }
+  if (num_vertices(G))
+    cout << "Error in the computation of feedback vertex set\n";
+  return reordered_vertices;
+}
+
+pair<int, vector<int>>
+VariableDependencyGraph::sortedStronglyConnectedComponents() const
+{
+  vector<int> vertex2scc(num_vertices(*this));
+  auto v_index = get(vertex_index, *this);
+
+  // Compute SCCs and create mapping from vertices to unordered SCCs
+  int num_scc = strong_components(static_cast<base>(*this), make_iterator_property_map(vertex2scc.begin(), v_index));
+
+  // Create directed acyclic graph (DAG) associated to the SCCs
+  adjacency_list<vecS, vecS, directedS> dag(num_scc);
+  for (int i = 0; i < static_cast<int>(num_vertices(*this)); i++)
+    {
+      out_edge_iterator it_out, out_end;
+      auto vi = vertex(i, *this);
+      for (tie(it_out, out_end) = out_edges(vi, *this); it_out != out_end; ++it_out)
+        {
+          int t_b = vertex2scc[v_index[target(*it_out, *this)]];
+          int s_b = vertex2scc[v_index[source(*it_out, *this)]];
+          if (s_b != t_b)
+            add_edge(s_b, t_b, dag);
+        }
+    }
+
+  /* Compute topological sort of DAG (ordered list of unordered SCC)
+     Note: the order is reversed. */
+  vector<int> reverseOrdered2unordered;
+  topological_sort(dag, back_inserter(reverseOrdered2unordered));
+
+  // Construct mapping from unordered SCC to ordered SCC
+  vector<int> unordered2ordered(num_scc);
+  for (int j = 0; j < num_scc; j++)
+    unordered2ordered[reverseOrdered2unordered[num_scc-j-1]] = j;
+
+  // Update the mapping of vertices to (now sorted) SCCs
+  for (int i = 0; i < static_cast<int>(num_vertices(*this)); i++)
+    vertex2scc[i] = unordered2ordered[vertex2scc[i]];
+
+  return { num_scc, vertex2scc };
+}
diff --git a/src/VariableDependencyGraph.hh b/src/VariableDependencyGraph.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6974c0206bd13f24f44019ad8fd7655db37d35f3
--- /dev/null
+++ b/src/VariableDependencyGraph.hh
@@ -0,0 +1,90 @@
+/*
+ * Copyright © 2009-2020 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/>.
+ */
+
+#ifndef _VARIABLEDEPENDENCYGRAPH_HH
+#define _VARIABLEDEPENDENCYGRAPH_HH
+
+#include <map>
+#include <vector>
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wold-style-cast"
+#include <boost/graph/adjacency_list.hpp>
+#pragma GCC diagnostic pop
+
+using namespace std;
+
+using VertexProperty_t = boost::property<boost::vertex_index_t, int,
+                                           boost::property<boost::vertex_index1_t, int,
+                                                           boost::property<boost::vertex_degree_t, int,
+                                                                           boost::property<boost::vertex_in_degree_t, int,
+                                                                                           boost::property<boost::vertex_out_degree_t, int>>>>>;
+
+/* Class used to store a graph representing dependencies between variables.
+   Used in the block decomposition. */
+class VariableDependencyGraph
+  : public boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>
+{
+public:
+  using color_t = map<boost::graph_traits<VariableDependencyGraph>::vertex_descriptor, boost::default_color_type>;
+  using base = boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>;
+
+  VariableDependencyGraph(int n);
+  //! Extracts a subgraph
+  /*!
+    \param[in] select_index The vertex indices to select
+    \return The subgraph
+
+    The property vertex_index1 of the subgraph contains indices of the original
+    graph.
+  */
+  VariableDependencyGraph extractSubgraph(const set<int> &select_index) const;
+  //! Return the feedback set
+  pair<VariableDependencyGraph, set<int>> minimalSetOfFeedbackVertices() const;
+  //! Reorder the recursive variables
+  /*! They appear first in a quasi triangular form and they are followed by the feedback variables */
+  vector<int> reorderRecursiveVariables(const set<int> &feedback_vertices) const;
+  /* Computes the strongly connected components (SCCs) of the graph, and sort them
+     topologically.
+     Returns the number of SCCs, and a mapping of vertex indices to sorted SCC
+     indices. */
+  pair<int, vector<int>> sortedStronglyConnectedComponents() const;
+private:
+  // Remove a vertex (including all edges to and from it); takes a vertex descriptor
+  void suppress(vertex_descriptor vertex_to_eliminate);
+  // Remove a vertex (including all edges to and from it); takes a vertex index
+  void suppress(int vertex_num);
+  /* Remove a vertex, but keeping the paths that go through it (i.e. by adding
+     edges that directly connect vertices that would otherwise be connected
+     through the vertex to be removed) */
+  void eliminate(vertex_descriptor vertex_to_eliminate);
+  // Internal helper for hasCycle()
+  bool hasCycleDFS(vertex_descriptor u, color_t &color, vector<int> &circuit_stack) const;
+  // Determine whether the graph has a cycle
+  bool hasCycle() const;
+  // Print on stdout a description of the graph
+  void print() const;
+  bool vertexBelongsToAClique(vertex_descriptor vertex) const;
+  bool eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree();
+  bool eliminationOfVerticesBelongingToAClique();
+  // The suppressed vertices are stored in feedback set
+  bool suppressionOfVerticesWithLoop(set<int> &feed_back_vertices);
+};
+
+#endif // _VARIABLEDEPENDENCYGRAPH_HH