From ea184312309acf11800dcb13895397367a82114b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 30 Jul 2020 15:11:09 +0200
Subject: [PATCH] Various improvements to mjdgges MEX

---
 mex/sources/mjdgges/mjdgges.F08 | 58 ++++++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 19 deletions(-)

diff --git a/mex/sources/mjdgges/mjdgges.F08 b/mex/sources/mjdgges/mjdgges.F08
index 43996e503..75c3f658d 100644
--- a/mex/sources/mjdgges/mjdgges.F08
+++ b/mex/sources/mjdgges/mjdgges.F08
@@ -1,3 +1,23 @@
+! Wrapper around LAPACK’s dgges (generalized Schur decomposition) that gives a
+! better access to error conditions than does MATLAB’s qz.
+!
+! Syntax:
+!   [ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
+!
+! Inputs:
+!   e            [double]  real square (n×n) matrix
+!   d            [double]  real square (n×n) matrix
+!   qz_criterium [double]  scalar (of the form 1+ε)
+!   zhreshold    [double]  used for detecting eigenvalues too close to 0÷0
+!
+! Outputs:
+!   ss           [double]  (n×n) quasi-triangular matrix
+!   tt           [double]  (n×n) quasi-triangular matrix
+!   zz           [double]  (n×n) orthogonal matrix
+!   sdim         [integer] scalar, number of stable eigenvalues
+!   eigval       [complex] (n×1) vector of generalized eigenvalues
+!   info         [integer] scalar, error code of dgges (or 30 if eigenvalue close to 0÷0)
+
 ! Copyright © 2006-2020 Dynare Team
 !
 ! This file is part of Dynare.
@@ -44,7 +64,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   type(c_ptr), dimension(*), intent(out) :: plhs
   integer(c_int), intent(in), value :: nlhs, nrhs
 
-  integer(c_size_t) :: m1, n1, m2, n2
+  integer(c_size_t) :: n
   real(real64) :: zhreshold
   integer(blint) :: n_bl, lwork, info_bl, sdim_bl
   real(real64), dimension(:), allocatable :: alpha_r, alpha_i, beta, work
@@ -61,41 +81,41 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
      return
   end if
 
-  m1 = mxGetM(prhs(1))
-  n1 = mxGetN(prhs(1))
-  m2 = mxGetM(prhs(2))
-  n2 = mxGetN(prhs(2))
-
+  n = mxGetM(prhs(1))
   if (.not. mxIsDouble(prhs(1)) .or. mxIsComplex(prhs(1)) &
       .or. .not. mxIsDouble(prhs(2)) .or. mxIsComplex(prhs(2)) &
-      .or. m1 /= n1 .or. m2 /= n1 .or. m2 /= n2) then
-     call mexErrMsgTxt("MJDGGES requires two square real matrices of the same dimension.")
+      .or. mxGetN(prhs(1)) /= n .or. mxGetM(prhs(2)) /= n .or. mxGetN(prhs(2)) /= n) then
+     call mexErrMsgTxt("MJDGGES: first two arguments should be real matrices of the same dimension")
      return
   end if
 
   ! Set criterium for stable eigenvalues
   if (nrhs >= 3 .and. mxGetM(prhs(3)) > 0) then
-     associate (crit_arg => mxGetPr(prhs(3)))
-       criterium = crit_arg(1)
-     end associate
+     if (.not. (mxIsScalar(prhs(3)) .and. mxIsNumeric(prhs(3)))) then
+        call mexErrMsgTxt("MJDGGES: third argument (qz_criterium) should be a numeric scalar")
+        return
+     end if
+     criterium = mxGetScalar(prhs(3))
   else
      criterium = 1_real64 + 1e-6_real64
   end if
 
   ! set criterium for 0/0 generalized eigenvalues */
   if (nrhs == 4 .and. mxGetM(prhs(4)) > 0) then
-     associate (zhresh_arg => mxGetPr(prhs(4)))
-       zhreshold = zhresh_arg(1)
-     end associate
+     if (.not. (mxIsScalar(prhs(4)) .and. mxIsNumeric(prhs(4)))) then
+        call mexErrMsgTxt("MJDGGES: fourth argument (zhreshold) should be a numeric scalar")
+        return
+     end if
+     zhreshold = mxGetScalar(prhs(4))
   else
      zhreshold = 1e-6_real64
   end if
 
-  plhs(1) = mxCreateDoubleMatrix(n1, n1, mxREAL)
-  plhs(2) = mxCreateDoubleMatrix(n1, n1, mxREAL)
-  plhs(3) = mxCreateDoubleMatrix(n1, n1, mxREAL)
+  plhs(1) = mxCreateDoubleMatrix(n, n, mxREAL)
+  plhs(2) = mxCreateDoubleMatrix(n, n, mxREAL)
+  plhs(3) = mxCreateDoubleMatrix(n, n, mxREAL)
   plhs(4) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
-  plhs(5) = mxCreateDoubleMatrix(n1, 1_mwSize, mxCOMPLEX)
+  plhs(5) = mxCreateDoubleMatrix(n, 1_mwSize, mxCOMPLEX)
   plhs(6) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
 
   s => mxGetPr(plhs(1))
@@ -117,7 +137,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
     t = b
   end associate
 
-  n_bl = int(n1, blint)
+  n_bl = int(n, blint)
   lwork = 16*n_bl + 16
   allocate(alpha_r(n_bl), alpha_i(n_bl), beta(n_bl), bwork(n_bl), work(lwork))