From e8700bc046d5fcea47ca8e42fb3782c53f9f62bd Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Thu, 9 Jul 2009 16:35:07 +0000
Subject: [PATCH] Implemented MATLAB part of the block_mfs option to steady.
 (Still need to fix a bug)

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2829 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/block_mfs_steadystate.m | 30 ++++++++++++++++++++++++++++++
 matlab/global_initialization.m |  6 +++++-
 matlab/homotopy1.m             | 10 +++-------
 matlab/homotopy2.m             |  9 +++------
 matlab/homotopy3.m             | 22 ++++++++++------------
 matlab/steady_.m               | 14 ++++++++++++++
 preprocessor/StaticModel.cc    | 15 +++++++--------
 7 files changed, 72 insertions(+), 34 deletions(-)
 create mode 100644 matlab/block_mfs_steadystate.m

diff --git a/matlab/block_mfs_steadystate.m b/matlab/block_mfs_steadystate.m
new file mode 100644
index 0000000000..f1c7de4df4
--- /dev/null
+++ b/matlab/block_mfs_steadystate.m
@@ -0,0 +1,30 @@
+function [r, g1] = block_mfs_steadystate(y, b)
+% Wrapper around the *_static.m file, for use with dynare_solve,
+% when block_mfs option is given to steady.
+
+% Copyright (C) 2009 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/>.
+
+  global M_ oo_
+  
+  ss = oo_.steady_state;
+  indx = M_.blocksMFS{b};
+  
+  ss(indx) = y;
+  x = [oo_.exo_steady_state; oo_.exo_det_steady_state];
+
+  eval(['[r,g1] = ' M_.fname '_static(b, ss, x, M_.params);']);
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 80107e87ab..8e3128b5fb 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -220,4 +220,8 @@ function global_initialization()
 
   % prior analysis
   options_.prior_mc = 20000;
-  options_.prior_analysis_endo_var_list = [];
\ No newline at end of file
+  options_.prior_analysis_endo_var_list = [];
+  
+  % block decomposition + minimum feedback set for steady state computation
+  options_.block_mfs = 0;
+  
\ No newline at end of file
diff --git a/matlab/homotopy1.m b/matlab/homotopy1.m
index 4c7241f340..1d76eda5d3 100644
--- a/matlab/homotopy1.m
+++ b/matlab/homotopy1.m
@@ -77,13 +77,9 @@ function homotopy1(values, step_nbr)
     oo_.exo_steady_state(values(ix,2)) = points(ix,i);
     oo_.exo_det_steady_state(values(ixd,2)) = points(ixd,i);
 
-    [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
-                                            oo_.steady_state,...
-                                            options_.jacobian_flag, ...	    
-                                            [oo_.exo_steady_state; ...
-                        oo_.exo_det_steady_state], M_.params);
-  
-    if check
+    try
+      steady_;
+    catch
       error('HOMOTOPY mode 1: failed')
     end
   end
diff --git a/matlab/homotopy2.m b/matlab/homotopy2.m
index 29078bcd52..2425e12c5a 100644
--- a/matlab/homotopy2.m
+++ b/matlab/homotopy2.m
@@ -99,13 +99,10 @@ function homotopy2(values, step_nbr)
       end
 
       disp([ 'HOMOTOPY mode 2: lauching solver with ' strtrim(varname) ' = ' num2str(v) ' ...'])
-      [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
-                                              oo_.steady_state,...
-                                              options_.jacobian_flag, ...	    
-                                              [oo_.exo_steady_state; ...
-                          oo_.exo_det_steady_state], M_.params);
       
-      if check
+      try
+        steady_;
+      catch
         error('HOMOTOPY mode 2: failed')
       end
     end
diff --git a/matlab/homotopy3.m b/matlab/homotopy3.m
index 1ee7dbedd3..c3fc9ddbe8 100644
--- a/matlab/homotopy3.m
+++ b/matlab/homotopy3.m
@@ -89,17 +89,10 @@ function homotopy3(values, step_nbr)
     oo_.exo_det_steady_state(values(ixd,2)) = curvalues(ixd);
     
     old_ss = oo_.steady_state;
-    [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
-                                            oo_.steady_state,...
-                                            options_.jacobian_flag, ...	    
-                                            [oo_.exo_steady_state; ...
-                        oo_.exo_det_steady_state], M_.params);
-  
-    if check
-      disp('HOMOTOPY mode 3: failed step, now dividing increment by 2...')
-      inc = inc/2;
-      oo_.steady_state = old_ss;
-    else
+
+    try
+      steady_;
+      
       if length([kplus; kminus]) == nv
         return
       end
@@ -110,7 +103,12 @@ function homotopy3(values, step_nbr)
       end
       oldvalues = curvalues;
       inc = 2*inc;
-    end
+    catch
+      disp('HOMOTOPY mode 3: failed step, now dividing increment by 2...')
+      inc = inc/2;
+      oo_.steady_state = old_ss;
+    end      
+      
     curvalues = oldvalues + inc;
     kplus = find(curvalues(iplus) >= targetvalues(iplus));
     curvalues(iplus(kplus)) = targetvalues(iplus(kplus));
diff --git a/matlab/steady_.m b/matlab/steady_.m
index 6c74f6ca22..a7c7c6a62d 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -62,6 +62,20 @@ function steady_()
 					      options_.jacobian_flag, ...	    
 					      [oo_.exo_steady_state;oo_.exo_det_steady_state],indv);
     end
+  elseif options_.block_mfs
+    for b = 1:size(M_.blocksMFS,1)
+      n = size(M_.blocksMFS{b}, 1);
+      ss = oo_.steady_state;
+      if n ~= 0
+        [y, check] = dynare_solve('block_mfs_steadystate', ...
+                                  ss(M_.blocksMFS{b}), ...
+                                  options_.jacobian_flag, b);
+        ss(M_.blocksMFS{b}) = y;
+      end
+      [r, g1, oo_.steady_state] = feval([M_.fname '_static'], b, ss, ...
+                          [oo_.exo_steady_state; ...
+                          oo_.exo_det_steady_state], M_.params);
+    end
   else
     [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
 				     oo_.steady_state,...
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index f36565d4b3..fd6f4d3334 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -517,15 +517,11 @@ StaticModel::writeOutput(ostream &output) const
   if (!block_mfs)
     return;
 
-  output << "M_.blocks = cell(" << blocks.size() << ", 1);" << endl
-         << "M_.blocksMFS = cell(" << blocksMFS.size() << ", 1);" << endl;
+  output << "M_.blocksMFS = cell(" << blocksMFS.size() << ", 1);" << endl;
   for(int b = 0; b < (int) blocks.size(); b++)
     {
-      output << "M_.blocks{" << b+1 << "} = [ ";
-      transform(blocks[b].begin(), blocks[b].end(), ostream_iterator<int>(output, " "), bind2nd(plus<int>(), 1));
-      output << "];" << endl
-             << "M_.blocksMFS{" << b+1 << "} = [ ";
-      transform(blocksMFS[b].begin(), blocksMFS[b].end(), ostream_iterator<int>(output, " "), bind2nd(plus<int>(), 1));
+      output << "M_.blocksMFS{" << b+1 << "} = [ ";
+      transform(blocksMFS[b].begin(), blocksMFS[b].end(), ostream_iterator<int>(output, "; "), bind2nd(plus<int>(), 1));
       output << "];" << endl;
     }
 }
@@ -533,7 +529,7 @@ StaticModel::writeOutput(ostream &output) const
 void
 StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) const
 {
-  output << "function [y, residual, g1] = " << func_name << "(nblock, y, x, params)" << endl
+  output << "function [residual, g1, y] = " << func_name << "(nblock, y, x, params)" << endl
          << "  switch nblock" << endl;
 
   for(int b = 0; b < (int) blocks.size(); b++)
@@ -596,4 +592,7 @@ StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) c
           i++;
         }
     }
+
+  output << "  end" << endl
+         << "end" << endl;
 }
-- 
GitLab