Commit ea184312 by Sébastien Villemot

### Various improvements to mjdgges MEX

parent f4a31a0d
 ! 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)) ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!