From c25ff09307dfae6bf158d2d1e77369eb761eebd9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 16 Jun 2020 12:43:56 +0200
Subject: [PATCH] =?UTF-8?q?Block=20decomposition:=20for=20=E2=80=9Csolve?=
 =?UTF-8?q?=20two=20boundaries=E2=80=9D=20block,=20no=20longer=20compute?=
 =?UTF-8?q?=20=E2=80=9Cb=E2=80=9D=20(-residuals+g1*y)=20in=20the=20dynamic?=
 =?UTF-8?q?=20function?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 matlab/lnsrch1_wrapper_two_boundaries.m |  2 +-
 matlab/solve_two_boundaries.m           | 14 +++++---------
 preprocessor                            |  2 +-
 3 files changed, 7 insertions(+), 11 deletions(-)

diff --git a/matlab/lnsrch1_wrapper_two_boundaries.m b/matlab/lnsrch1_wrapper_two_boundaries.m
index 0da783520..d1c5ecfcc 100644
--- a/matlab/lnsrch1_wrapper_two_boundaries.m
+++ b/matlab/lnsrch1_wrapper_two_boundaries.m
@@ -47,5 +47,5 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 
 %reshape the input arguments of the dynamic function
 y(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
-[r, T, g1, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, y_size);
+[r, T, g1]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, y_size);
 ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*y_size, 1);
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index e6d5c77d9..f638975ad 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -78,16 +78,14 @@ ilu_setup.milu = 'off';
 ilu_setup.thresh = 1;
 ilu_setup.udiag = 0;
 max_resa=1e100;
-Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
-g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
 reduced = 0;
 while ~(cvg==1 || iter>maxit_)
-    [r, T, g1, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, Blck_size);
+    [r, T, g1]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, Blck_size);
     preconditioner = 2;
     g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
-    term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
-    term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
-    b = b - term1 - term2;
+    ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';
+    ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
+    b=-ra+g1a*ya;
     [max_res, max_indx]=max(max(abs(r')));
     if ~isreal(r)
         max_res = (-max_res^2)^0.5;
@@ -120,7 +118,7 @@ while ~(cvg==1 || iter>maxit_)
                                 disp(['    trying to correct the Jacobian matrix:']);
                                 disp(['    correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
                             end
-                            dx = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya;
+                            dx = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya_save;
                             y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
                             continue
                         else
@@ -158,7 +156,6 @@ while ~(cvg==1 || iter>maxit_)
                 end
             end
         end
-        ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';
         ya_save=ya;
         g1aa=g1a;
         ba=b;
@@ -299,7 +296,6 @@ while ~(cvg==1 || iter>maxit_)
                 end
             end
         elseif stack_solve_algo==4
-            ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
             stpmx = 100 ;
             stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);
             nn=1:size(ra,1);
diff --git a/preprocessor b/preprocessor
index 849697937..85b900363 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 849697937579eb432cbf507e1b6e5f0493e26091
+Subproject commit 85b900363d58b42aa0d0dc4d5e3a2d02bf792448