From 5670d3938fd89200b66449ce255516be62fcaa3c Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Mon, 27 Jun 2022 15:29:48 +0200
Subject: [PATCH] dynare_solve: deal with Dulmage-Mendelsohn decomposition
returns a non-square block
Related to https://git.dynare.org/Dynare/dynare/-/issues/1851
(manually cherry picked from commit 5788f1bc7146b5d7be654c2ddefe5e59cb5158a2)
---
matlab/dynare_solve.m | 28 +++++++++++++++++++++++-----
1 file changed, 23 insertions(+), 5 deletions(-)
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 31e427d78b..3e323b2c92 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -311,11 +311,29 @@ elseif ismember(options.solve_algo, [2, 12, 4])
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
end
end
- [x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
- options.gstep, ...
- tolf, options.solve_tolx, ...
- maxit, options.debug, args{:});
- fre = true;
+ blockcolumns=s(i+1)-s(i);
+ if blockcolumns ~= blocklength
+ %non-square-block in DM; check whether initial value is solution
+ [fval_check, fjac] = feval(f, x, args{:});
+ if norm(fval_check(j1(j))) < tolf
+ errorflag = false;
+ errorcode = 0;
+ continue
+ end
+ end
+ if blockcolumns>=blocklength
+ %(under-)determined block
+ [x, errorflag, errorcode] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
+ options.gstep, ...
+ tolf, options.solve_tolx, maxit, ...
+ options.debug, args{:});
+ fre = true;
+ else
+ fprintf('\nDYNARE_SOLVE (solve_algo=2|4|12): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.\n')
+ %overdetermined block
+ errorflag = true;
+ errorcode = 0;
+ end
if errorflag
return
end
--
GitLab