m4/ax_mexopts.m4

@@ -33,7 +33,7 @@ case ${MATLAB_ARCH} in
     MATLAB_DEFS="-D_GNU_SOURCE -DNDEBUG"
     MATLAB_CFLAGS="-fexceptions -fPIC -pthread"
     MATLAB_CXXFLAGS="-fPIC -pthread"
-    MATLAB_FCFLAGS="-fPIC -fexceptions"
+    MATLAB_FCFLAGS="-fPIC -fexceptions -g"
     MATLAB_LDFLAGS_NOMAP="-shared -Wl,--no-undefined -Wl,-rpath-link,$MATLAB/bin/glnxa64 -L$MATLAB/bin/glnxa64"
     MATLAB_LDFLAGS="$MATLAB_LDFLAGS_NOMAP -Wl,--version-script,$MATLAB/extern/lib/glnxa64/mexFunction.map"
     MATLAB_LIBS="-lmx -lmex -lmat -lm -lmwlapack -lmwblas"

mex/build/k_order_perturbation_fortran.am

+mex_PROGRAMS = k_order_perturbation_fortran 
+
+k_order_perturbation_fortran_FCFLAGS = $(AM_FCFLAGS) -I../libkordersim -L/home/normann/Dropbox/dynare/SuiteSparse/lib
+
+nodist_k_order_perturbation_fortran_SOURCES = \
+	mexFunction.F08
+
+k_order_perturbation_fortran_LDADD = ../libkordersim/libkordersim.a -lcholmod -lspqr
+
+BUILT_SOURCES = $(nodist_k_order_perturbation_fortran_SOURCES)
+CLEANFILES = $(nodist_k_order_perturbation_fortran_SOURCES)
+
+%.F08: $(top_srcdir)/../../sources/k_order_perturbation_fortran/%.F08
+	$(LN_S) -f $< $@

mex/build/libkordersim.am

@@ -8,7 +8,8 @@ nodist_libkordersim_a_SOURCES = \
 	partitions.f08 \
 	simulation.f08 \
 	struct.f08 \
-	pthread.F08
+	pthread.F08 \
+	cholmod.F08
 
 
 BUILT_SOURCES = $(nodist_libkordersim_a_SOURCES)

mex/build/matlab/Makefile.am

 ACLOCAL_AMFLAGS = -I ../../../m4
 
-SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update
+SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim k_order_perturbation_fortran local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update
 
 # libdynare++ must come before gensylv and k_order_perturbation
 if ENABLE_MEX_DYNAREPLUSPLUS

mex/build/matlab/configure.ac

@@ -160,6 +160,7 @@ AC_CONFIG_FILES([Makefile
 	             sobol/Makefile
 		         local_state_space_iterations/Makefile
                  libkordersim/Makefile
+                 k_order_perturbation_fortran/Makefile
                  folded_to_unfolded_dr/Makefile
                  k_order_simul/Makefile
                  k_order_mean/Makefile

mex/build/matlab/k_order_perturbation_fortran/Makefile.am

+include ../mex.am
+include ../../k_order_perturbation_fortran.am

mex/build/matlab/mex.am

@@ -49,3 +49,8 @@ pthread.mod: pthread.o
 
 pthread.F08: $(top_srcdir)/../../sources/pthread.F08
 	$(LN_S) -f $< $@
+
+cholmod.mod: cholmod.o
+
+cholmod.F08: $(top_srcdir)/../../sources/cholmod.F08
+	$(LN_S) -f $< $@

mex/build/octave/configure.ac

@@ -163,6 +163,7 @@ AC_CONFIG_FILES([Makefile
 		         sobol/Makefile
 		         local_state_space_iterations/Makefile
                  libkordersim/Makefile 
+                 k_order_perturbation_fortran/Makefile
                  folded_to_unfolded_dr/Makefile
                  k_order_simul/Makefile
                  k_order_mean/Makefile

mex/build/octave/k_order_perturbation_fortran/Makefile.am

+EXEEXT = .mex
+include ../mex.am
+include ../../k_order_perturbation_fortran.am

mex/build/octave/mex.am

@@ -57,3 +57,8 @@ pthread.mod: pthread.o
 
 pthread.F08: $(top_srcdir)/../../sources/pthread.F08
 	$(LN_S) -f $< $@
+
+cholmod.mod: cholmod.o
+
+cholmod.F08: $(top_srcdir)/../../sources/cholmod.F08
+	$(LN_S) -f $< $@

mex/sources/Makefile.am

@@ -7,6 +7,7 @@ EXTRA_DIST = \
 	defines.F08 \
 	matlab_mex.F08 \
 	pthread.F08 \
+	cholmod.F08 \
 	mjdgges \
 	kronecker \
 	bytecode \
@@ -17,6 +18,7 @@ EXTRA_DIST = \
 	sobol \
 	local_state_space_iterations \
 	libkordersim \
+	k_order_perturbation_fortran \
 	folded_to_unfolded_dr \
 	k_order_simul \
    k_order_mean \

mex/sources/cholmod.F08

+! Fortran interface for a subset of CHOLMOD functions
+
+! Copyright © 2023 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+module cholmod
+   use iso_fortran_env
+   use iso_c_binding
+   implicit none
+   
+
+   integer(c_int), parameter :: CHOLMOD_REAL = 1
+   integer(c_int), parameter :: CHOLMOD_ZOMPLEX = 3
+
+   integer(c_int), parameter :: CHOLMOD_INT = 0
+
+   integer(c_int), parameter :: CHOLMOD_DOUBLE = 0
+   integer(c_int), parameter :: CHOLMOD_LONG = 2
+
+   integer(c_int), parameter :: TRUE = 1
+   integer(c_int), parameter :: FALSE = 0
+
+   integer(c_int), parameter :: SPQR_ORDERING_FIXED = 0
+
+   real(c_double), parameter :: SPQR_DEFAULT_TOL = -2.
+
+   type, bind(c) :: cholmod_sparse_struct
+      integer(c_size_t) :: nrow, ncol, nzmax
+      type(c_ptr) :: p, i, nz, x, z
+      integer(c_int) :: stype, itype, xtype, dtype, sorted, packed
+   end type cholmod_sparse_struct
+
+   type, bind(c) :: cholmod_triplet_struct
+      integer(c_size_t) :: nrow, ncol, nzmax, nnz
+      type(c_ptr) :: i, j, x, z
+      integer(c_int) :: stype, itype, xtype, dtype
+   end type cholmod_triplet_struct
+
+   type, bind(c) :: cholmod_method_struct
+      real(c_double) :: lnz, fl, prune_dense, prune_dense2, nd_oksep, &
+                       &other_1(4)
+      integer(c_size_t) :: nd_small, other_2(4) 
+      integer(c_int) :: aggressive, order_for_lu, nd_compress, nd_camd, &
+                       &nd_components, ordering
+      integer(c_size_t) :: other_3(4) 
+   end type cholmod_method_struct
+
+   type, bind(c) :: cholmod_common_struct
+      real(c_double) :: dbound, grow0, grow1
+      integer(c_size_t) :: grow2, maxrank
+      real(c_double) :: supernodal_switch
+      integer(c_int) :: supernodal, final_asis, final_super, final_ll, &
+                       &final_pack, final_monotonic, final_resymbol
+      real(c_double) :: zrelax(3)
+      integer(c_size_t) :: nrelax(3)
+      integer(c_int) :: prefer_zomplex, prefer_upper, &
+                       &quick_return_if_not_posdef, prefer_binary, &
+                       &print, precise, try_catch
+      type(c_funptr) :: error_handler
+      integer(c_int) :: nmethods, current, selected
+      type(cholmod_method_struct) :: method(10)
+      integer(c_int) :: postorder, default_nesdis
+      real(c_double) :: metis_memory, metis_dswitch
+      integer(c_size_t) :: metis_nswitch, nrow
+      integer(c_int64_t) :: mark
+      integer(c_size_t) :: iworksize, xworksize
+      type(c_ptr) :: Flag, Head, Xwork, Iwork
+      integer(c_int) :: itype, dtype, no_workspace_reallocate, status
+      real(c_double) :: fl, lnz, anz, modfl
+      integer(c_size_t) :: malloc_count, memory_usage, memory_inuse
+      real(c_double) :: nrealloc_col, nrealloc_factor, ndbounds_hit, &
+                       &rowfacfl, aatfl
+      integer(c_int) :: called_nd, blas_ok
+      real(c_double) :: SPQR_grain, SPQR_small
+      integer(c_int) :: SPQR_shrink, SPQR_nthreads
+      real(c_double) :: SPQR_flopcount, SPQR_analyze_time, &
+                       &SPQR_factorize_time, SPQR_solve_time, &
+                       &SPQR_flopcount_bound, SPQR_tol_used, &
+                       &SPQR_norm_E_fro
+      integer(c_int64_t) :: SPQR_istat(10)
+      integer(c_int) :: useGPU
+      integer(c_size_t) :: maxGpuMemBytes
+      real(c_double) :: maxGpuMemFraction
+      integer(c_size_t) :: gpuMemorySize
+      real(c_double) :: gpuKernelTime
+      integer(c_int64_t) :: gpuFlops
+      integer(c_int) :: gpuNumKernelLaunches
+      type(c_ptr) :: cublasHandle, gpuStream(8), &
+                    &cublasEventPotrf (3), updateCKernelsComplete, &
+                    &updateCBuffersFree(8), &
+                    &dev_mempool
+      integer(c_size_t) :: dev_mempool_size
+      type(c_ptr) :: host_pinned_mempool
+      integer(c_size_t) :: host_pinned_mempool_size, devBuffSize
+      integer(c_int) :: ibuffer
+      real(c_double) :: syrkStart, cholmod_cpu_gemm_time, &
+                       &cholmod_cpu_syrk_time, cholmod_cpu_trsm_time, &
+                       &cholmod_cpu_potrf_time, cholmod_gpu_gemm_time, &
+                       &cholmod_gpu_syrk_time, cholmod_gpu_trsm_time, &
+                       &cholmod_gpu_potrf_time, cholmod_assemble_time, &
+                       &cholmod_assemble_time2
+      integer(c_size_t) :: cholmod_cpu_gemm_calls, cholmod_cpu_syrk_calls, &
+                          &cholmod_cpu_trsm_calls, cholmod_cpu_potrf_calls, &
+                          &cholmod_gpu_gemm_calls, cholmod_gpu_syrk_calls, &
+                          &cholmod_gpu_trsm_calls, cholmod_gpu_potrf_calls
+      real(c_double) :: chunk
+      integer(c_int) :: nthreads_max
+   end type cholmod_common_struct
+
+   interface
+      integer(c_int) function c_cholmod_start(common) bind(c, name="cholmod_start")
+         use iso_c_binding
+         implicit none 
+         type(c_ptr), value :: common
+      end function c_cholmod_start
+
+      integer(c_int) function c_cholmod_finish(common) bind(c, name="cholmod_finish")
+         use iso_c_binding
+         implicit none
+         type(c_ptr), value :: common
+      end function c_cholmod_finish
+
+      type(c_ptr) function c_cholmod_submatrix(A, rset, rsize, cset, csize, &
+                                              &values, sorted, common) &
+                                              &bind(c, name="cholmod_submatrix")
+         use iso_c_binding
+         implicit none
+         type(c_ptr), value, intent(in) :: A, rset, cset, common
+         integer(c_int64_t), value, intent(in) :: rsize, csize
+         integer(c_int), value, intent(in) :: values, sorted
+      end function c_cholmod_submatrix
+
+      integer(c_int64_t) function c_qr(order, tol, econ, A, Q, R, E, common) bind(c, name="SuiteSparseQR_C_QR")
+         use iso_c_binding
+         implicit none      
+         integer(c_int), value, intent(in) :: order
+         real(c_double), value, intent(in) :: tol
+         integer(c_int64_t), value, intent(in) :: econ
+         type(c_ptr), intent(in) :: Q, R
+         type(c_ptr), value, intent(in) :: A, E, common
+      end function c_qr
+   end interface
+
+contains
+
+   integer(c_int) function cholmod_start(cholmod_common)
+      type(cholmod_common_struct), target :: cholmod_common
+      cholmod_start = c_cholmod_start(c_loc(cholmod_common))
+   end function cholmod_start
+
+   integer(c_int) function cholmod_finish(cholmod_common)
+      type(cholmod_common_struct), target :: cholmod_common
+      cholmod_finish = c_cholmod_finish(c_loc(cholmod_common))
+   end function cholmod_finish
+
+   ! cholmod_submatrix:  C = A (r,c), where r and c are arbitrary vectors
+   ! A         matrix to subreference
+   ! rset      set of row indices, duplicates OK
+   ! rsize     size of r; rsize < 0 denotes ":"
+   ! cset	   set of column indices, duplicates OK
+   ! csize 	   size of c; csize < 0 denotes ":"
+   ! values	   if TRUE compute the numerical values of C */
+   ! sorted    if TRUE then return C with sorted columns */
+   ! rsize < 0 denotes ":" in MATLAB notation, or more precisely 0:(A->nrow)-1.
+   ! In this case, r can be NULL.  An rsize of zero, or r = NULL and rsize >= 0,
+   ! denotes "[ ]" in MATLAB notation (the empty set).
+   ! Similar rules hold for csize.
+   function cholmod_submatrix_col(A, cset, common)
+      type(cholmod_sparse_struct), target, intent(in) :: A
+      integer(c_int32_t), dimension(:), target, intent(in) :: cset
+      type(cholmod_common_struct), target, intent(inout) :: common
+      type(cholmod_sparse_struct), pointer :: cholmod_submatrix_col 
+      call c_f_pointer(c_cholmod_submatrix(c_loc(A), c_null_ptr, -1_c_int64_t, &
+                                          &c_loc(cset),                        &
+                                          &int(size(cset), c_int64_t),         &
+                                          &TRUE, TRUE,                       &
+                                          &c_loc(common)),                     &
+                                          &cholmod_submatrix_col)
+   end function cholmod_submatrix_col
+
+   integer(c_int64_t) function qr(A, Q, R, common)
+      type(cholmod_sparse_struct), target, intent(inout) :: A, Q, R
+      type(cholmod_common_struct), target, intent(inout) :: common
+      qr = c_qr(SPQR_ORDERING_FIXED, SPQR_DEFAULT_TOL,      &
+               &int(A%nrow, c_int64_t), c_loc(A), c_loc(Q), &
+               &c_loc(R), c_null_ptr, c_loc(common))
+   end function qr
+
+! int64_t SuiteSparseQR_C_QR          // returns rank(A) estimate, (-1) if failure
+! (
+!     // inputs:
+!     int ordering,           // all, except 3:given treated as 0:fixed
+!     double tol,             // columns with 2-norm <= tol are treated as 0
+!     int64_t econ,              // e = max(min(m,econ),rank(A))
+!     cholmod_sparse *A,      // m-by-n sparse matrix to factorize
+!     // outputs:
+!     cholmod_sparse **Q,     // m-by-e sparse matrix
+!     cholmod_sparse **R,     // e-by-n sparse matrix
+!     int64_t **E,               // size n column permutation, NULL if identity
+!     cholmod_common *cc      // workspace and parameters
+! )
+
+
+end module cholmod

mex/sources/k_order_perturbation_fortran/mexFunction.F08

+! Copyright © 2022 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+!
+!  input:
+!       oo_
+!       M_
+!       options_
+module usefulfun
+   use iso_c_binding
+   use matlab_mex
+   use cholmod
+   implicit none
+
+   integer(c_int), parameter :: CHOLMOD_OK = 0
+
+contains
+
+   subroutine c_chars_to_f_string(c_string, f_string)
+      character(len=1,kind=c_char), dimension(*), intent(in) :: c_string
+      character(len=*), intent(out) :: f_string
+      integer :: i
+      i=1
+      do while ((c_string(i) /= c_null_char) .and. (i <= len(f_string)))
+         f_string(i:i) = c_string(i)
+         i=i+1
+      end do
+      if (i<len(f_string)) f_string(i:) = ' '
+   end subroutine c_chars_to_f_string
+
+   subroutine error_matlab(status, file, line, message) bind(c)
+      integer(c_int), value, intent(in) :: status, line
+      character(kind=c_char), dimension(*), intent(in) :: file, message
+      character(kind=c_char, len=100) :: line_str, status_str, file_str, message_str
+      call c_chars_to_f_string(message, message_str)
+      call c_chars_to_f_string(file, file_str)
+      write (line_str, "(i0)") line
+      write (status_str, "(i0)") status
+      if (status < CHOLMOD_OK) then
+         call mexPrintf ("ERROR: file "//trim(file_str)//" line "//     &
+                        &trim(line_str)//", status "//trim(status_str)//&
+                        &new_line(status_str))
+      end if
+      call mexErrMsgTxt(trim(message_str))
+   end subroutine error_matlab
+
+   subroutine fill_sparse_matlab(A, A_mx)
+      type(cholmod_sparse_struct), intent(inout) :: A
+      type(c_ptr), intent(in) :: A_mx
+      integer(c_size_t), pointer, contiguous :: jc(:), ir(:)
+      integer(c_int32_t), dimension(:), allocatable, target :: jc32, ir32
+      integer :: i
+      if (.not. (c_associated(A_mx) .and. mxIsSparse(A_mx) &
+         &.and. (.not. mxIsComplex(A_mx)))) then
+         call mexErrMsgTxt("The input should be a real sparse MATLAB matrix.")
+      end if
+      A%nrow = mxGetM(A_mx)
+      A%ncol = mxGetN(A_mx)
+      jc => mxGetJc(A_mx)
+      jc32 = [ (int(jc(i), c_int32_t), i=1, size(jc)) ]
+      A%p = c_loc(jc32)
+      ir => mxGetIr(A_mx)
+      ir32 = [ (int(ir(i), c_int32_t), i=1, size(ir)) ]
+      A%i = c_loc(ir32)
+      A%nzmax = int(jc(A%ncol+1), c_size_t)
+      A%packed = TRUE
+      A%sorted = TRUE
+      A%nz = c_null_ptr
+      A%itype = CHOLMOD_LONG
+      A%dtype = CHOLMOD_DOUBLE
+      A%stype = 0_c_int
+#if MX_HAS_INTERLEAVED_COMPLEX
+      A%x = mxGetDoubles_internal(A_mx)
+#else
+      A%x = mxGetPr_internal(A_mx)
+#endif
+      A%xtype = CHOLMOD_REAL
+   end subroutine fill_sparse_matlab
+
+end module usefulfun
+
+subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
+   use iso_fortran_env
+   use iso_c_binding
+   use struct
+   use cholmod
+   use matlab_mex
+   use usefulfun
+   implicit none
+
+   type(c_ptr), dimension(*), intent(in), target :: prhs
+   type(c_ptr), dimension(*), intent(out) :: plhs
+   integer(c_int), intent(in), value :: nlhs, nrhs
+   type(c_ptr) :: oo_mx, dr_mx, M_mx, options_mx, fname_mx, &  
+                 &steady_state_mx, params_mx, y_mx, x_mx, &
+                 &rowval_mx, colval_mx, colptr_mx, ys_mx, &
+                 &g1_mx
+   character(kind=c_char, len=:), allocatable :: fname
+   character(kind=c_char, len=100) :: funcname
+   integer :: endo_nbr, exo_nbr, order, i, retval
+   real(real64), pointer, contiguous :: ys(:), params(:)
+   real(real64), allocatable :: ext_steady_state(:), x(:)
+   type(c_ptr), allocatable :: call_rhs(:), call_lhs(:)
+   type(cholmod_common_struct), target :: cholmod_common
+   integer(int32), pointer, contiguous :: rowval(:), colval(:)
+   integer(c_int) :: retval_cholmod
+   type(cholmod_sparse_struct) :: g1, Q, R
+   type(cholmod_sparse_struct), pointer :: H_theta
+   integer(c_int64_t) :: rank
+   integer(c_int32_t), dimension(:), allocatable :: indcs
+   integer(c_size_t), dimension(:), pointer, contiguous :: jc
+
+   oo_mx = prhs(1)
+   M_mx = prhs(2)
+   options_mx = prhs(3)
+
+   ! 0. Check the input arguments and import the necessary elements
+   ! Check the number of input arguments
+   if (nrhs /= 3) then
+      call mexErrMsgTxt("Must have 3 input parameters")
+   end if
+
+   if (.not. mxIsStruct(oo_mx)) then
+      call mexErrMsgTxt("1st argument (oo) should be a struct")
+   end if
+   if (.not. mxIsStruct(M_mx)) then
+      call mexErrMsgTxt("2nd argument (M) should be a struct")
+   end if
+   if (.not. mxIsStruct(options_mx)) then
+      call mexErrMsgTxt("3rd argument (options) should be a struct")
+   end if
+
+   ! Import the necessary elements from dr, M and options
+   order = get_int_field(options_mx, "order")
+   endo_nbr = get_int_field(M_mx, "endo_nbr")
+   exo_nbr = get_int_field(M_mx, "exo_nbr")
+
+   dr_mx = mxGetField(oo_mx, 1_mwIndex, "dr")
+   if (.not. (c_associated(dr_mx) .and. mxIsStruct(dr_mx))) then
+      call mexErrMsgTxt("oo_.dr should be a struct")
+   end if
+
+   fname_mx = mxGetField(M_mx, 1_mwIndex, "fname")
+   if (.not. (c_associated(fname_mx) .and. mxIsChar(fname_mx) &
+       .and. mxGetM(fname_mx) == 1)) then
+      call mexErrMsgTxt("M_.fname should be a character string")
+   end if
+   fname = mxArrayToString(fname_mx)
+
+   params_mx = mxGetField(M_mx, 1_mwIndex, "params")
+   if (.not. (c_associated(params_mx) .and. mxIsDouble(params_mx) .and. &
+       (.not. mxIsComplex(params_mx)) .and. (.not. mxIsSparse(params_mx)))) then
+      call mexErrMsgTxt("M_.params should be a real dense array")
+   end if
+
+   ys_mx = mxGetField(dr_mx, 1_mwIndex, "ys")
+   if (.not. (c_associated(ys_mx) .and. mxIsDouble(ys_mx) .and. &
+       (.not. mxIsComplex(ys_mx)) .and. (.not. mxIsSparse(ys_mx)))) then
+      call mexErrMsgTxt("dr.ys should be a real dense array")
+   end if
+#if MX_HAS_INTERLEAVED_COMPLEX
+   ys => mxGetDoubles(ys_mx)
+#else
+   ys => mxGetPr(ys_mx)
+#endif
+
+   steady_state_mx = mxGetField(oo_mx, 1_mwIndex, "steady_state")
+   if (.not. (c_associated(steady_state_mx) .and. mxIsDouble(steady_state_mx) &
+      &.and. (.not. mxIsComplex(steady_state_mx)) .and. &
+       (.not. mxIsSparse(steady_state_mx)))) then
+      call mexErrMsgTxt("oo_.steady_state should be a real dense array")
+   end if
+
+   rowval_mx = mxGetField(M_mx, 1_mwIndex, "dynamic_g1_sparse_rowval")
+   if (.not. (c_associated(rowval_mx) .and. mxIsInt32(rowval_mx) .and. &
+       (.not. mxIsSparse(rowval_mx)))) then
+      call mexErrMsgTxt("M_.colval should be an integer dense vector")
+   end if
+
+   colval_mx = mxGetField(M_mx, 1_mwIndex, "dynamic_g1_sparse_colval")
+   if (.not. (c_associated(colval_mx) .and. mxIsInt32(colval_mx) .and. &
+       (.not. mxIsSparse(colval_mx)))) then 
+      call mexErrMsgTxt("M_.colval should be an integer dense vector")
+   end if
+
+   colptr_mx = mxGetField(M_mx, 1_mwIndex, "dynamic_g1_sparse_colptr")
+   if (.not. (c_associated(colptr_mx) .and. mxIsInt32(colptr_mx) .and. &
+       (.not. mxIsSparse(colptr_mx)))) then 
+      call mexErrMsgTxt("M_.colptr should be an integer dense vector")
+   end if
+
+! #if MX_HAS_INTERLEAVED_COMPLEX
+!    rowval => mxGetInt32s(rowval_mx)
+!    colval => mxGetInt32s(colval_mx)
+! #else
+!    rowval => mxGetData(rowval_mx)
+!    colval => mxGetData(colval_mx)
+! #endif
+
+   ! 1. Call to the Matlab <model>_dynamic routine to get model derivatives
+   ! (i) Building the proper input variables
+   ! Vector of endogenous variables at times t-1, t and t+1 denoted y
+   ! The endogenous variables take the value of their deterministic 
+   ! steady-state value stored in oo_.dr.ys
+   allocate(ext_steady_state(3*endo_nbr))
+   ext_steady_state(1:endo_nbr) = ys
+   ext_steady_state(endo_nbr+1:2*endo_nbr) = ys
+   ext_steady_state(2*endo_nbr+1:3*endo_nbr) = ys
+   y_mx = mxCreateDoubleMatrix(int(size(ext_steady_state), mwSize), 1_mwSize, mxREAL)
+   mxGetPr(y_mx) = ext_steady_state
+
+   ! Vector of exogenous variables, which are set to zero
+   allocate(x(exo_nbr))
+   x = 0.
+   x_mx = mxCreateDoubleMatrix(int(exo_nbr,mwSize), 1_mwSize, mxREAL)
+   mxGetPr(x_mx) = x
+
+   ! The vector of steady-state values as eventually specified in 
+   ! the STEADY_STATE block and stored in oo_.steady_state
+   ! is kept in steady_state_mx. Similarly, the other inputs are already 
+   ! imported.
+   
+   ! (ii) Calling the Matlab <model>_dynamic routine and storing the
+   ! relevant outputs in Fortran variables
+   allocate(call_lhs(3))
+   call_rhs = [y_mx, x_mx, params_mx, steady_state_mx, rowval_mx, colval_mx, colptr_mx]
+   write (funcname, "(a,'.sparse.dynamic_g', i0)")  fname, 1
+   retval = mexCallMATLAB(3_c_int, call_lhs, 7_c_int, call_rhs, trim(funcname))
+   if (retval /= 0) then
+      call mexErrMsgTxt("Trouble calling "//trim(funcname))
+   end if
+
+   g1_mx = call_lhs(1)
+   if (.not. (c_associated(g1_mx) .and. mxIsSparse(g1_mx))) then
+      call mexErrMsgTxt(trim(funcname)//" 1st output variable should be a sparse array")
+   end if
+
+   ! do i=1,3
+   !    call mxDestroyArray(call_lhs(i)[ (int(i, c_int32_t), i=1*endo_nbr-1, 1*endo_nbr-1) ])
+   ! end do
+   ! deallocate(call_lhs)
+
+   ! 2. Sparse routines
+   ! Initializing the CHOLMOD workspace
+   retval_cholmod = cholmod_start(cholmod_common)
+   if (retval_cholmod /= 1_c_int) then
+      call mexErrMsgTxt("The initialization of the CHOLMOD workspace failed.")
+   end if
+   cholmod_common%error_handler = c_funloc(error_matlab)
+
+   ! Anderson, G.
+   ! "Solving Linear Rational Expectations Models: A Horse Race"
+   ! Computational Economics, 2008, vol. 31, issue 2, pages 95-113
+   !
+   ! Anderson, G.
+   ! "A Reliable and Computationally Efficient Algorithm for Imposing the
+   ! Saddle Point Property in Dynamic Models"
+   ! Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
+   ! pages 472-489
+   call fill_sparse_matlab(g1, g1_mx)  
+   ! (i) State-space transition matrix and auxiliary initial conditions
+   jc => mxGetJc(g1_mx)
+   
+   ! Extract H_θ
+   indcs = [ (int(i, c_int32_t), i=2*endo_nbr, 3*endo_nbr-1) ]
+   H_theta => cholmod_submatrix_col(g1, indcs, cholmod_common)
+
+   rank = qr(H_theta, Q, R, cholmod_common)
+
+   retval_cholmod = cholmod_finish(cholmod_common)
+   if (retval_cholmod /= 1_c_int) then
+      call mexErrMsgTxt("The deallocation of the CHOLMOD workspace failed.")
+   end if
+
+contains
+
+
+end subroutine mexFunction
\ No newline at end of file

mex/sources/matlab_mex.F08

@@ -153,6 +153,11 @@ module matlab_mat
        type(c_ptr), intent(in), value :: pm
      end function mxGetPr_internal
      
+     logical(c_bool) function mxIsInt32(pm) bind(c, name="mxIsInt32"//API_VER)
+        use iso_c_binding
+        type(c_ptr), intent(in), value :: pm
+     end function mxIsInt32
+
      ! Noncomplex integer
 #if MX_HAS_INTERLEAVED_COMPLEX
      type(c_ptr) function mxGetInt32s_internal(pa) bind(c, name="mxGetInt32s"//API_VER)
@@ -187,15 +192,15 @@ module matlab_mat
        type(c_ptr), intent(in), value :: pm
      end function mxIsSparse
 
-     type(c_ptr) function mxGetIr(pm) bind(c, name="mxGetIr"//API_VER2)
+     type(c_ptr) function mxGetIr_internal(pm) bind(c, name="mxGetIr"//API_VER2)
        use iso_c_binding
        type(c_ptr), intent(in), value :: pm
-     end function mxGetIr
+     end function mxGetIr_internal
 
-     type(c_ptr) function mxGetJc(pm) bind(c, name="mxGetJc"//API_VER2)
+     type(c_ptr) function mxGetJc_internal(pm) bind(c, name="mxGetJc"//API_VER2)
        use iso_c_binding
        type(c_ptr), intent(in), value :: pm
-     end function mxGetJc
+     end function mxGetJc_internal
 
      ! Nonnumeric types
      type(c_ptr) function mxGetData(pm) bind(c, name="mxGetData"//API_VER)
@@ -345,6 +350,18 @@ contains
   end function mxGetPi
 #endif
 
+  function mxGetIr(pm)
+    type(c_ptr), intent(in) :: pm
+    integer(mwIndex), dimension(:), pointer, contiguous :: mxGetIr
+    call c_f_pointer(mxGetIr_internal(pm), mxGetIr, [ mxGetNumberOfElements(pm) ])
+  end function mxGetIr
+
+  function mxGetJc(pm)
+    type(c_ptr), intent(in) :: pm
+    integer(mwIndex), dimension(:), pointer, contiguous :: mxGetJc
+    call c_f_pointer(mxGetJc_internal(pm), mxGetJc, [ int(mxGetN(pm)+1) ])
+  end function mxGetJc
+
   function mxGetLogicals(array_ptr)
     type(c_ptr), intent(in) :: array_ptr
     logical(mxLogical), dimension(:), pointer, contiguous :: mxGetLogicals

preprocessor @ 389a2647

-Subproject commit 114d8eadfb3dee24ab65f291253076c4d0ea83c3
+Subproject commit 389a2647d3b1067d3af54596ad3bd076fab8a4b4