From 16eabbbc4e9691518f4bf46d68d75fa35cc2a6ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Fri, 22 Jul 2022 12:40:41 +0200 Subject: [PATCH] block_trust_region MEX: improve treatment of non-square blocks in Dulmage-Mendelsohn decomposition MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit – before erroring out, check whether the residuals for the block are already zero (in which case, move to next block) – improve error message that is printed otherwise Note that trying to solve under-determined blocks (as in dynare_solve.m) would require too many changes in the existing code, so let’s leave it out. Closes: #1851 --- mex/sources/block_trust_region/mexFunction.f08 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08 index 1092d4233..59da65547 100644 --- a/mex/sources/block_trust_region/mexFunction.f08 +++ b/mex/sources/block_trust_region/mexFunction.f08 @@ -201,14 +201,24 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') end if end if + if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then + ! Non-square block in DM decomposition + ! Before erroring out, check whether we are not already at the solution for this block + ! See also #1851 + if (norm2(fvec(blocks(i)%row_indices)) < tolf) then + cycle + else + call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13|14): 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.") + end if + end if + block real(real64), dimension(size(blocks(i)%col_indices)) :: x_block x_indices => blocks(i)%col_indices f_indices => blocks(i)%row_indices x_all => x - if (size(x_indices) /= size(f_indices)) then - call mexErrMsgTxt("Non-square block") - end if x_block = x(x_indices) call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor) x(x_indices) = x_block