From 10fdc42516f4982cba30f8bc6ac159814a1d48f7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 30 Nov 2022 14:47:48 +0100
Subject: [PATCH] block_trust_region MEX: add support for sparse Jacobian

If the function to be solved returns a sparse Jacobian, simply convert it to a
dense representation.
---
 .../block_trust_region/matlab_fcn_closure.F08      | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/mex/sources/block_trust_region/matlab_fcn_closure.F08 b/mex/sources/block_trust_region/matlab_fcn_closure.F08
index 403214fd54..398ac0b417 100644
--- a/mex/sources/block_trust_region/matlab_fcn_closure.F08
+++ b/mex/sources/block_trust_region/matlab_fcn_closure.F08
@@ -134,11 +134,21 @@ contains
 
     ! Handle Jacobian
     if (present(fjac)) then
-       if (.not. mxIsDouble(call_lhs(2)) .or. mxIsSparse(call_lhs(2)) &
+       if (.not. mxIsDouble(call_lhs(2)) &
             .or. mxGetM(call_lhs(2)) /= n_all .or. mxGetN(call_lhs(2)) /= n_all) &
-            call mexErrMsgTxt("Second output argument of the function must be a dense double float&
+            call mexErrMsgTxt("Second output argument of the function must be a double float&
             & square matrix whose dimension matches the length of the first input argument")
 
+       ! If Jacobian is sparse, convert it to a dense matrix
+       if (mxIsSparse(call_lhs(2))) then
+          block
+            type(c_ptr) :: dense_jacobian
+            if (mexCallMATLAB(1_c_int, dense_jacobian, 1_c_int, call_lhs(2), "full") /= 0) &
+                 call mexErrMsgTxt("Error converting Jacobian from sparse to dense representation")
+            call_lhs(2) = dense_jacobian
+          end block
+       end if
+
        if (.not. mxIsComplex(call_lhs(2))) then ! Real case
           block
             real(real64), dimension(:,:), pointer, contiguous :: fjac_all
-- 
GitLab