From cab65dabb57cf6f0ede1dad3fb5df9d95f129c76 Mon Sep 17 00:00:00 2001 From: NormannR Date: Fri, 27 Aug 2021 13:55:51 +0200 Subject: [PATCH] Implements `local_state_space_iteration_fortran`, the Fortran replacement of `local_state_space_iteration_k` (Ref #1802) - Creates the library `libkordersim` with all the relevant Fortran routines to `folded_to_unfolded_dr` and `local_state_space_iteration_fortran` - Implements `folded_to_unfolded_dr`, which converts folded decision rule matrices to their unfolded counterparts --- mex/build/folded_to_unfolded_dr.am | 14 + mex/build/libkordersim.am | 29 ++ .../local_state_space_iteration_fortran.am | 14 + mex/build/matlab/Makefile.am | 2 +- mex/build/matlab/configure.ac | 3 + .../matlab/folded_to_unfolded_dr/Makefile.am | 2 + mex/build/matlab/libkordersim/Makefile.am | 2 + .../Makefile.am | 2 + mex/build/octave/Makefile.am | 2 +- mex/build/octave/configure.ac | 3 + .../octave/folded_to_unfolded_dr/Makefile.am | 3 + mex/build/octave/libkordersim/Makefile.am | 3 + .../Makefile.am | 3 + mex/sources/Makefile.am | 3 + .../folded_to_unfolded_dr/mexFunction.f08 | 101 +++++++ mex/sources/libkordersim/partitions.f08 | 256 ++++++++++++++++++ mex/sources/libkordersim/pascal.f08 | 89 ++++++ mex/sources/libkordersim/simulation.f08 | 127 +++++++++ mex/sources/libkordersim/sort.f08 | 39 +++ mex/sources/libkordersim/struct.f08 | 36 +++ .../mexFunction.f08 | 173 ++++++++++++ mex/sources/matlab_mex.F08 | 17 ++ .../local_state_space_iteration_k_test.mod | 36 +++ 23 files changed, 957 insertions(+), 2 deletions(-) create mode 100644 mex/build/folded_to_unfolded_dr.am create mode 100644 mex/build/libkordersim.am create mode 100644 mex/build/local_state_space_iteration_fortran.am create mode 100644 mex/build/matlab/folded_to_unfolded_dr/Makefile.am create mode 100644 mex/build/matlab/libkordersim/Makefile.am create mode 100644 mex/build/matlab/local_state_space_iteration_fortran/Makefile.am create mode 100644 mex/build/octave/folded_to_unfolded_dr/Makefile.am create mode 100644 mex/build/octave/libkordersim/Makefile.am create mode 100644 mex/build/octave/local_state_space_iteration_fortran/Makefile.am create mode 100644 mex/sources/folded_to_unfolded_dr/mexFunction.f08 create mode 100644 mex/sources/libkordersim/partitions.f08 create mode 100644 mex/sources/libkordersim/pascal.f08 create mode 100644 mex/sources/libkordersim/simulation.f08 create mode 100644 mex/sources/libkordersim/sort.f08 create mode 100644 mex/sources/libkordersim/struct.f08 create mode 100644 mex/sources/local_state_space_iteration_fortran/mexFunction.f08 diff --git a/mex/build/folded_to_unfolded_dr.am b/mex/build/folded_to_unfolded_dr.am new file mode 100644 index 000000000..c96e8a009 --- /dev/null +++ b/mex/build/folded_to_unfolded_dr.am @@ -0,0 +1,14 @@ +mex_PROGRAMS = folded_to_unfolded_dr + +folded_to_unfolded_dr_FCFLAGS = $(AM_FCFLAGS) -Warray-temporaries -I../libkordersim + +nodist_folded_to_unfolded_dr_SOURCES = \ + mexFunction.f08 + +folded_to_unfolded_dr_LDADD = ../libkordersim/libkordersim.a + +BUILT_SOURCES = $(nodist_folded_to_unfolded_dr_SOURCES) +CLEANFILES = $(nodist_folded_to_unfolded_dr_SOURCES) + +%.f08: $(top_srcdir)/../../sources/folded_to_unfolded_dr/%.f08 + $(LN_S) -f $< $@ diff --git a/mex/build/libkordersim.am b/mex/build/libkordersim.am new file mode 100644 index 000000000..192bf3292 --- /dev/null +++ b/mex/build/libkordersim.am @@ -0,0 +1,29 @@ +noinst_LIBRARIES = libkordersim.a + +nodist_libkordersim_a_SOURCES = \ + blas_lapack.F08 \ + matlab_mex.F08 \ + pascal.f08 \ + sort.f08 \ + partitions.f08 \ + simulation.f08 \ + struct.f08 + +BUILT_SOURCES = $(nodist_libkordersim_a_SOURCES) +CLEANFILES = $(nodist_libkordersim_a_SOURCES) + +struct.o: matlab_mex.mod +struct.mod: struct.o + +pascal.mod: pascal.o + +sort.mod: sort.o + +partitions.o: sort.mod pascal.mod +partitions.mod: partitions.o + +simulation.o: partitions.mod lapack.mod +simulation.mod: simulation.o + +%.f08: $(top_srcdir)/../../sources/libkordersim/%.f08 + $(LN_S) -f $< $@ diff --git a/mex/build/local_state_space_iteration_fortran.am b/mex/build/local_state_space_iteration_fortran.am new file mode 100644 index 000000000..e3ee6eee7 --- /dev/null +++ b/mex/build/local_state_space_iteration_fortran.am @@ -0,0 +1,14 @@ +mex_PROGRAMS = local_state_space_iteration_fortran + +local_state_space_iteration_fortran_FCFLAGS = $(AM_FCFLAGS) -Warray-temporaries -I../libkordersim + +nodist_local_state_space_iteration_fortran_SOURCES = \ + mexFunction.f08 + +local_state_space_iteration_fortran_LDADD = ../libkordersim/libkordersim.a + +BUILT_SOURCES = $(nodist_local_state_space_iteration_fortran_SOURCES) +CLEANFILES = $(nodist_local_state_space_iteration_fortran_SOURCES) + +%.f08: $(top_srcdir)/../../sources/local_state_space_iteration_fortran/%.f08 + $(LN_S) -f $< $@ diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am index 58f376cdc..ae0ca758b 100644 --- a/mex/build/matlab/Makefile.am +++ b/mex/build/matlab/Makefile.am @@ -4,7 +4,7 @@ SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_ if ENABLE_MEX_DYNAREPLUSPLUS -SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations +SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran endif if ENABLE_MEX_MS_SBVAR diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac index 69468b0d5..cfec90594 100644 --- a/mex/build/matlab/configure.ac +++ b/mex/build/matlab/configure.ac @@ -162,6 +162,9 @@ AC_CONFIG_FILES([Makefile block_kalman_filter/Makefile sobol/Makefile local_state_space_iterations/Makefile + libkordersim/Makefile + folded_to_unfolded_dr/Makefile + local_state_space_iteration_fortran/Makefile perfect_foresight_problem/Makefile num_procs/Makefile block_trust_region/Makefile diff --git a/mex/build/matlab/folded_to_unfolded_dr/Makefile.am b/mex/build/matlab/folded_to_unfolded_dr/Makefile.am new file mode 100644 index 000000000..b6ea3f920 --- /dev/null +++ b/mex/build/matlab/folded_to_unfolded_dr/Makefile.am @@ -0,0 +1,2 @@ +include ../mex.am +include ../../folded_to_unfolded_dr.am \ No newline at end of file diff --git a/mex/build/matlab/libkordersim/Makefile.am b/mex/build/matlab/libkordersim/Makefile.am new file mode 100644 index 000000000..727b6b607 --- /dev/null +++ b/mex/build/matlab/libkordersim/Makefile.am @@ -0,0 +1,2 @@ +include ../mex.am +include ../../libkordersim.am diff --git a/mex/build/matlab/local_state_space_iteration_fortran/Makefile.am b/mex/build/matlab/local_state_space_iteration_fortran/Makefile.am new file mode 100644 index 000000000..448c76eb2 --- /dev/null +++ b/mex/build/matlab/local_state_space_iteration_fortran/Makefile.am @@ -0,0 +1,2 @@ +include ../mex.am +include ../../local_state_space_iteration_fortran.am \ No newline at end of file diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index a669cf839..93295839b 100644 --- a/mex/build/octave/Makefile.am +++ b/mex/build/octave/Makefile.am @@ -4,7 +4,7 @@ SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_ if ENABLE_MEX_DYNAREPLUSPLUS -SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations +SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran endif if ENABLE_MEX_MS_SBVAR diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index 6a38305d0..9d361b95b 100644 --- a/mex/build/octave/configure.ac +++ b/mex/build/octave/configure.ac @@ -144,6 +144,9 @@ AC_CONFIG_FILES([Makefile block_kalman_filter/Makefile sobol/Makefile local_state_space_iterations/Makefile + libkordersim/Makefile + folded_to_unfolded_dr/Makefile + local_state_space_iteration_fortran/Makefile perfect_foresight_problem/Makefile num_procs/Makefile block_trust_region/Makefile diff --git a/mex/build/octave/folded_to_unfolded_dr/Makefile.am b/mex/build/octave/folded_to_unfolded_dr/Makefile.am new file mode 100644 index 000000000..d0cad16e3 --- /dev/null +++ b/mex/build/octave/folded_to_unfolded_dr/Makefile.am @@ -0,0 +1,3 @@ +EXEEXT = .mex +include ../mex.am +include ../../folded_to_unfolded_dr.am \ No newline at end of file diff --git a/mex/build/octave/libkordersim/Makefile.am b/mex/build/octave/libkordersim/Makefile.am new file mode 100644 index 000000000..2488da73b --- /dev/null +++ b/mex/build/octave/libkordersim/Makefile.am @@ -0,0 +1,3 @@ +EXEEXT = .mex +include ../mex.am +include ../../libkordersim.am diff --git a/mex/build/octave/local_state_space_iteration_fortran/Makefile.am b/mex/build/octave/local_state_space_iteration_fortran/Makefile.am new file mode 100644 index 000000000..8dcce8dc6 --- /dev/null +++ b/mex/build/octave/local_state_space_iteration_fortran/Makefile.am @@ -0,0 +1,3 @@ +EXEEXT = .mex +include ../mex.am +include ../../local_state_space_iteration_fortran.am diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am index 1840bbf6a..4a7ce297f 100644 --- a/mex/sources/Makefile.am +++ b/mex/sources/Makefile.am @@ -16,6 +16,9 @@ EXTRA_DIST = \ block_kalman_filter \ sobol \ local_state_space_iterations \ + libkordersim \ + folded_to_unfolded_dr \ + local_state_space_iteration_fortran \ gensylv \ dynare_simul_ \ perfect_foresight_problem \ diff --git a/mex/sources/folded_to_unfolded_dr/mexFunction.f08 b/mex/sources/folded_to_unfolded_dr/mexFunction.f08 new file mode 100644 index 000000000..7e772abd4 --- /dev/null +++ b/mex/sources/folded_to_unfolded_dr/mexFunction.f08 @@ -0,0 +1,101 @@ +! Copyright © 2021 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 . + +subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') + use iso_fortran_env + use iso_c_binding + use struct + use simulation + use matlab_mex + use partitions + 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) :: M_mx, options_mx, dr_mx, tmp, g + type(pol), dimension(:), allocatable, target :: fdr + integer :: order, npred, nboth, nstatic, nfwrd, endo_nbr, exo_nbr, nys, nvar + type(pascal_triangle) :: p + type(uf_matching), dimension(:), allocatable :: matching + character(kind=c_char, len=10), dimension(:), allocatable :: fieldnames + integer :: d, m, n + + dr_mx = prhs(1) + M_mx = prhs(2) + options_mx = prhs(3) + + ! Checking the consistence and validity of input arguments + if (nrhs /= 3 .or. nlhs /= 1) then + call mexErrMsgTxt("Must have exactly 3 inputs and 1 output") + end if + if (.not. mxIsStruct(dr_mx)) then + call mexErrMsgTxt("1st argument (dr) 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 + + nstatic = get_int_field(M_mx, "nstatic") + npred = get_int_field(M_mx, "npred") + nboth = get_int_field(M_mx, "nboth") + nfwrd = get_int_field(M_mx, "nfwrd") + endo_nbr = nstatic+npred+nboth+nfwrd + exo_nbr = get_int_field(M_mx, "exo_nbr") + order = get_int_field(options_mx, "order") + nys = npred+nboth + nvar = nys+exo_nbr + + allocate(fdr(0:order), fieldnames(0:order)) + do d = 0, order + write (fieldnames(d), '(a2, i1)') "g_", d + tmp = mxGetField(dr_mx, 1_mwIndex, trim(fieldnames(d))) + if (.not. (c_associated(tmp) .and. mxIsDouble(tmp))) then + call mexErrMsgTxt(trim(fieldnames(d))//" is not allocated in dr") + end if + m = int(mxGetM(tmp)) + n = int(mxGetN(tmp)) + allocate(fdr(d)%g(m,n)) + fdr(d)%g = reshape(mxGetPr(tmp), [m,n]) + end do + + plhs(1) = mxCreateStructMatrix(1_mwSize, 1_mwSize, fieldnames) + g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), 1_mwSize, mxREAL) + mxGetPr(g) = reshape(fdr(0)%g, [size(fdr(0)%g)]) + call mxSetField(plhs(1), 1_mwIndex, "g_0", g) + g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nvar, mwSize), mxREAL) + mxGetPr(g) = reshape(fdr(1)%g, [size(fdr(1)%g)]) + call mxSetField(plhs(1), 1_mwIndex, "g_1", g) + + if (order > 1) then + ! Compute the useful binomial coefficients from Pascal's triangle + p = pascal_triangle(nvar+order-1) + allocate(matching(2:order)) + ! Pinpointing the corresponding offsets between folded and unfolded tensors + do d=2,order + allocate(matching(d)%folded(nvar**d)) + call fill_folded_indices(matching(d)%folded, nvar, d, p) + g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nvar**d, mwSize), mxREAL) + mxGetPr(g) = reshape(fdr(d)%g(:,matching(d)%folded), [size(fdr(d)%g(:,matching(d)%folded))]) + call mxSetField(plhs(1), 1_mwIndex, trim(fieldnames(d)), g) + end do + end if + +end subroutine mexFunction \ No newline at end of file diff --git a/mex/sources/libkordersim/partitions.f08 b/mex/sources/libkordersim/partitions.f08 new file mode 100644 index 000000000..0fb0829e4 --- /dev/null +++ b/mex/sources/libkordersim/partitions.f08 @@ -0,0 +1,256 @@ +! Provides subroutines to manipulate indexes representing elements of +! a partition for a given integer +! i.e. elements p = (α₁,…,αₘ) where each αᵢ ∈ { 0, ..., n-1 } +! +! Copyright © 2021 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 . + +module partitions + use pascal + use sort + implicit none + + ! index represents the aforementioned (α₁,…,αₘ) objects + type index + integer, dimension(:), allocatable :: ind + end type index + + interface index + module procedure :: init_index + end interface index + + ! a dictionary that matches folded indices with folded offsets + type dict + integer :: pr ! pointer to the last added element in indices and offsets + type(index), dimension(:), allocatable :: indices ! list of folded indices + integer, dimension(:), allocatable :: offsets ! list of the associated offsets + end type dict + + interface dict + module procedure :: init_dict + end interface dict + + interface operator(/=) + module procedure :: diff_indices + end interface operator(/=) + + ! A type to contain the correspondence unfolded and folded offsets i.e. + ! folded(i) shall contain the folded offset corresponding to the unfolded offset i + type uf_matching + type(integer), dimension(:), allocatable :: folded + end type + +contains + + ! Constructor for the index type + type(index) function init_index(d, ind) + integer, intent(in) :: d + integer, dimension(d), intent(in) :: ind + allocate(init_index%ind(d)) + init_index%ind = ind + end function init_index + + ! Comparison for the index type. Returns true if the two indices are different + type(logical) function diff_indices(i1,i2) + type(index), intent(in) :: i1, i2 + if (size(i1%ind) /= size(i2%ind) .or. any(i1%ind /= i2%ind)) then + diff_indices = .true. + else + diff_indices = .false. + end if + end function diff_indices + + ! Constructor of the dict type + type(dict) function init_dict(n, d, p) + integer, intent(in) :: n, d + type(pascal_triangle), intent(in) :: p + integer :: size + size = get(d, n+d-1, p) + allocate(init_dict%indices(size), init_dict%offsets(size)) + init_dict%pr = 0 + end function init_dict + + ! Count the number of coordinates similar to the the first one for a given index + type(integer) function get_prefix_length(idx, d) + integer, intent(in) :: d + type(index), intent(in) :: idx + integer :: i + i = 1 + if (d>1) then + do while ((i < d) .and. (idx%ind(i+1) == idx%ind(1))) + i = i+1 + end do + end if + get_prefix_length = i + end function get_prefix_length + + ! Gets the folded index associated with an unfolded index + type(index) function u_index_to_f_index(idx, d) + type(index), intent(in) :: idx + integer, intent(in) :: d + u_index_to_f_index = index(d, idx%ind) + call sort_int(u_index_to_f_index%ind) + end function u_index_to_f_index + + ! Converts the offset of an unfolded tensor to the associated unfolded tensor index + ! Note that the index (α₁,…,αₘ) is such that αᵢ ∈ { 0, ..., n-1 } + ! and the offset is such that j ∈ {1, ..., nᵈ} + type(index) function u_offset_to_u_index(j, n, d) + integer, intent(in) :: j, n, d ! offset, number of variables and dimensions respectively + integer :: i, tmp, r + allocate(u_offset_to_u_index%ind(d)) + tmp = j-1 ! We substract 1 as j ∈ {1, ..., n} so that tmp ∈ {0, ..., n-1} and our modular operations work + do i=d,1,-1 + r = mod(tmp, n) + u_offset_to_u_index%ind(i) = r + tmp = (tmp-r)/n + end do + end function u_offset_to_u_index + + ! Converts a folded tensor index to the associated folded tensor offset + ! See the explanation in dynare++/tl/cc/tensor.cc for the function FTensor::getOffsetRecurse + ! Note that the index (α₁,…,αₘ) is such that αᵢ ∈ { 0, ..., n-1 } + ! and the offset is such that j ∈ {1, ..., ⎛n+d-1⎞ } + ! ⎝ d ⎠ + recursive function f_index_to_f_offset(idx, n, d, p) result(j) + type(index), intent(in) :: idx ! folded index + integer, intent(in) :: n, d ! number of variables and dimensions + type(pascal_triangle) :: p ! Pascal's triangle containing the relevant binomial coefficients + integer :: j, prefix + type(index) :: tmp + if (d == 0) then + j = 1 + else + prefix = get_prefix_length(idx,d) + tmp = index(d-prefix, idx%ind(prefix+1:) - idx%ind(1)) + j = get(d, n+d-1, p) - get(d, n-idx%ind(1)+d-1, p) + f_index_to_f_offset(tmp, n-idx%ind(1), d-prefix, p) + end if + end function f_index_to_f_offset + + ! Function that searches a value in an array of a given length + type(integer) function find(a, v, l) + integer, intent(in) :: l ! length of the array + type(index), dimension(l), intent(in) :: a ! array of indices + type(index) :: v ! element to be found + integer :: i + if (l == 0) then + find = 0 + else + i = 1 + do while (i <= l .and. a(i) /= v) + i = i+1 + end do + if (i == l+1) then + find = 0 + else + find = i + end if + end if + end function find + + ! Fills the folded offset array: + ! folded(i) shall contain the folded offset corresponding to the unfolded offset i + ! For each unfolded tensor offset + ! (a) compute the associated unfolded index (u_offset_to_u_index) + ! (b) compute the associated folded index (u_index_to_f_index) + ! (c) has the folded offset already been computed ? + ! (i) If yes, get the corresponding offset + ! (ii) If no, compute it (f_index_to_f_offset) and store it for reuse and as a result + subroutine fill_folded_indices(folded, n, d, p) + integer, intent(in) :: n, d + integer, dimension(n**d), intent(inout) :: folded + type(pascal_triangle), intent(in) :: p + type(dict) :: c + type(index) :: tmp + integer :: j, found + c = dict(n, d, p) + do j=1,n**d + tmp = u_offset_to_u_index(j,n,d) + tmp = u_index_to_f_index(tmp, d) + found = find(c%indices, tmp, c%pr) + if (found == 0) then + c%pr = c%pr+1 + c%indices(c%pr) = tmp + c%offsets(c%pr) = f_index_to_f_offset(tmp,n,d,p) + folded(j) = c%offsets(c%pr) + else + folded(j) = c%offsets(found) + end if + end do + end subroutine fill_folded_indices + + end module partitions + +! gfortran -o partitions partitions.f08 pascal.f08 sort.f08 +! ./partitions +! program test +! use partitions +! use pascal +! implicit none +! type(index) :: uidx, fidx, i1, i2 +! integer, dimension(:), allocatable :: folded +! integer :: i, uj, n, d +! type(pascal_triangle) :: p +! ! Unfolded indices and offsets +! ! 0,0,0 1 1,0,0 10 2,0,0 19 +! ! 0,0,1 2 1,0,1 11 2,0,1 20 +! ! 0,0,2 3 1,0,2 12 2,0,2 21 +! ! 0,1,0 4 1,1,0 13 2,1,0 22 +! ! 0,1,1 5 1,1,1 14 2,1,1 23 +! ! 0,1,2 6 1,1,2 15 2,1,2 24 +! ! 0,2,0 7 1,2,0 16 2,2,0 25 +! ! 0,2,1 8 1,2,1 17 2,2,1 26 +! ! 0,2,2 9 1,2,2 18 2,2,2 27 + +! ! Folded indices and offsets +! ! 0,0,0 1 1,1,1 7 2,2,2 10 +! ! 0,0,1 2 1,1,2 8 +! ! 0,0,2 3 1,2,2 9 +! ! 0,1,1 4 +! ! 0,1,2 5 +! ! 0,2,2 6 + +! n = 3 +! d = 3 +! uj = 8 +! p = pascal_triangle(n+d-1) + +! ! u_offset_to_u_index +! uidx = u_offset_to_u_index(uj,n,d) +! print '(3i2)', (uidx%ind(i), i=1,d) ! should display 0 2 1 + +! ! f_index_to_f_offset +! fidx = u_index_to_f_index(uidx, d) +! print '(i2)', f_index_to_f_offset(fidx, n, d, p) ! should display 5 + +! ! /= +! i1 = index(5, (/1,2,3,4,5/)) +! i2 = index(5, (/1,2,3,4,6/)) +! if (i1 /= i2) then +! print *, "Same!" +! else +! print *, "Different!" +! end if + +! ! fill_folded_indices +! allocate(folded(n**d)) +! call fill_folded_indices(folded,n,d) +! print *, "Matching offsets unfolded -> folded" +! print '(1000i4)', (i, i=1,n**d) +! print '(1000i4)', (folded(i), i=1,n**d) + +! end program test \ No newline at end of file diff --git a/mex/sources/libkordersim/pascal.f08 b/mex/sources/libkordersim/pascal.f08 new file mode 100644 index 000000000..bec24515e --- /dev/null +++ b/mex/sources/libkordersim/pascal.f08 @@ -0,0 +1,89 @@ +! Provides a subroutine to build Pascal's triangle +! +! Copyright © 2021 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 . + +module pascal + implicit none + + type line + integer, dimension(:), allocatable :: coeffs + end type line + + type pascal_triangle + integer :: d + type(line), dimension(:), allocatable :: lines + end type pascal_triangle + + interface pascal_triangle + module procedure :: init_pascal_triangle + end interface pascal_triangle + +contains + + ! Fills a pascal_triangle structure with the associated coefficients up to a given order + type(pascal_triangle) function init_pascal_triangle(d) + integer, intent(in) :: d + integer :: i, j + init_pascal_triangle%d = d + allocate(init_pascal_triangle%lines(d)) + ! Initializing the first line + allocate(init_pascal_triangle%lines(1)%coeffs(2)) + init_pascal_triangle%lines(1)%coeffs = 1 + ! Iterating Pascal's triangle formula + if (d > 1) then + do i=2,d + allocate(init_pascal_triangle%lines(i)%coeffs(i+1)) + init_pascal_triangle%lines(i)%coeffs(1) = 1 + init_pascal_triangle%lines(i)%coeffs(i+1) = 1 + do j=2,i + init_pascal_triangle%lines(i)%coeffs(j) = init_pascal_triangle%lines(i-1)%coeffs(j-1) & + + init_pascal_triangle%lines(i-1)%coeffs(j) + end do + end do + end if + end function init_pascal_triangle + + ! Returns ⎛n⎞ stored in pascal_triangle p + ! ⎝k⎠ + + type(integer) function get(k,n,p) + integer, intent(in) :: k, n + type(pascal_triangle), intent(in) :: p + get = p%lines(n)%coeffs(k+1) + end function get + +end module pascal + +! gfortran -o pascal pascal.f08 +! ./pascal +! program test +! use pascal +! type(pascal_triangle) :: p +! integer :: d +! read *, d +! p = pascal_triangle(d) +! do n=1,d +! do k=0,n +! if (k < n) then +! write (*,'(i2," ")', advance="no") get(k,n,p) +! else +! write (*,'(i2," ")') get(k,n,p) +! end if +! end do +! end do +! end program test \ No newline at end of file diff --git a/mex/sources/libkordersim/simulation.f08 b/mex/sources/libkordersim/simulation.f08 new file mode 100644 index 000000000..387f09bcf --- /dev/null +++ b/mex/sources/libkordersim/simulation.f08 @@ -0,0 +1,127 @@ +! Necessary routines and functions to carry out simulations +! +! A first step is to get the associated +! Copyright © 2021 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 . + +module simulation + use iso_fortran_env + use partitions + use lapack + implicit none + + ! Used to store the folded decision rule tensors + type :: pol + real(real64), allocatable :: g(:,:) + end type pol + + ! Used to store the different contracted tensors used in the Horner algorithm + ! type :: horner + ! real(real64), pointer, contiguous :: c(:,:), c_flat(:) + ! end type horner + type :: horner + real(real64), allocatable :: c(:,:) + end type horner + + +contains + + ! ! With MATMUL + ! ! Horner evaluation of the polynomial with coefficients stored in udr at the point dyu + ! subroutine eval(h, dyu, udr, ny, nvar, order) + ! type(horner), dimension(0:order), intent(inout) :: h + ! real(real64), dimension(nvar), intent(in) :: dyu + ! type(pol), intent(in) :: udr(0:order) + ! integer, intent(in) :: ny, nvar, order + ! integer :: d + ! if (order == 1) then + ! h(1)%c = udr(1)%g + ! else + ! call contract(h(order-1)%c, udr(order)%g, dyu, ny, nvar, order) + ! end if + ! do d=order-1,1,-1 + ! if (d > 1) then + ! h(d)%c = h(d)%c + udr(d)%g + ! call contract(h(d-1)%c, h(d)%c, dyu, ny, nvar, d) + ! else + ! h(d)%c = h(d)%c + udr(1)%g + ! end if + ! end do + ! h(0)%c(:,1) = matmul(h(1)%c, dyu) + udr(0)%g(:,1) + ! end subroutine eval + + ! ! Contracts the unfolded tensor t with respect to the vector x and stores the result in c + ! subroutine contract(c, t, x, nrows, nvar, d) + ! integer, intent(in) :: nrows, nvar, d + ! real(real64), dimension(nrows, nvar**d), intent(in) :: t + ! real(real64), dimension(nvar), intent(in) :: x + ! real(real64), dimension(nrows, nvar**(d-1)), intent(inout) :: c(:,:) + ! real(real64), dimension(nrows) :: tmp + ! integer :: i + ! do i=1,nvar**(d-1) + ! tmp = matmul(t(:,(i-1)*nvar+1:i*nvar), x) + ! c(:,i) = tmp + ! end do + ! end subroutine contract + + ! With DGEMV + ! Modifies y such that y := alpha*A*x + beta*y, where alpha, beta are scalars, + ! x and y are vectors, A is a m-by-n matrix + subroutine mult_vec(alpha, A, x, beta, y) + real(real64), intent(in) :: alpha, beta, x(:), A(:,:) + real(real64), intent(inout) :: y(:) + call dgemv("N", int(size(A,1), blint), int(size(A,2), blint), alpha, A, int(size(A,1), blint), x, 1_blint, beta, y, 1_blint) + end subroutine mult_vec + + ! Horner evaluation of the polynomial with coefficients stored in udr at the point dyu + subroutine eval(h, dyu, udr, ny, nvar, order) + type(horner), dimension(0:order), intent(inout) :: h + real(real64), dimension(nvar), intent(in) :: dyu + type(pol), intent(in) :: udr(0:order) + integer, intent(in) :: ny, nvar, order + integer :: d + if (order == 1) then + h(1)%c = udr(1)%g + else + call contract(h(order-1)%c, udr(order)%g, dyu, ny, nvar, order) + end if + do d=order-1,1,-1 + if (d > 1) then + h(d)%c = h(d)%c + udr(d)%g + call contract(h(d-1)%c, h(d)%c, dyu, ny, nvar, d) + else + h(d)%c = h(d)%c + udr(1)%g + end if + end do + h(0)%c = udr(0)%g + call mult_vec(1.0_real64, h(1)%c, dyu, 1.0_real64, h(0)%c(:,1)) + end subroutine eval + + ! Contracts the unfolded tensor t with respect to the vector x and stores the result in c + subroutine contract(c, t, x, nrows, nvar, d) + integer, intent(in) :: nrows, nvar, d + real(real64), dimension(nrows, nvar**d), intent(in) :: t + real(real64), dimension(nvar), intent(in) :: x + real(real64), dimension(nrows, nvar**(d-1)), intent(inout) :: c + integer :: i + do i=1,nvar**(d-1) + call mult_vec(1.0_real64, t(:,(i-1)*nvar+1:i*nvar), x, 0.0_real64, c(:,i)) + end do + end subroutine contract + +end module simulation + diff --git a/mex/sources/libkordersim/sort.f08 b/mex/sources/libkordersim/sort.f08 new file mode 100644 index 000000000..225f37852 --- /dev/null +++ b/mex/sources/libkordersim/sort.f08 @@ -0,0 +1,39 @@ +! Provides a subroutine to sort integer arrays in ascending order +! As the addressed arrays are small, I use the insertion sort algorithm +! +! Copyright © 2021 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 . + +module sort + implicit none + +contains + subroutine sort_int(l) + integer, dimension(:), intent(inout) :: l + integer :: i, j, x + do i=2,size(l) + x = l(i) + j = i + do while (j > 1 .and. l(j-1) > x) + l(j) = l(j-1) + j = j-1 + end do + l(j) = x + end do + end subroutine sort_int + +end module sort \ No newline at end of file diff --git a/mex/sources/libkordersim/struct.f08 b/mex/sources/libkordersim/struct.f08 new file mode 100644 index 000000000..c2a8e0337 --- /dev/null +++ b/mex/sources/libkordersim/struct.f08 @@ -0,0 +1,36 @@ +! Copyright © 2021 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 . +module struct + use iso_fortran_env + use iso_c_binding + use matlab_mex + implicit none + + contains + + type(integer) function get_int_field(struct, field) + type(c_ptr), intent(in) :: struct + character(*), intent(in) :: field + type(c_ptr) :: tmp + tmp = mxGetField(struct, 1_mwIndex, field) + if (.not. (c_associated(tmp) .and. mxIsScalar(tmp) .and. mxIsNumeric(tmp))) then + call mexErrMsgTxt("Field "//field//" should be a numeric scalar") + end if + get_int_field = int(mxGetScalar(tmp)) + end function get_int_field + +end module struct \ No newline at end of file diff --git a/mex/sources/local_state_space_iteration_fortran/mexFunction.f08 b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08 new file mode 100644 index 000000000..81195ec83 --- /dev/null +++ b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08 @@ -0,0 +1,173 @@ +! Copyright © 2021 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 . + +! input: +! order the order of approximation, needs order+1 derivatives +! nstat +! npred +! nboth +! nforw +! nexog +! ystart starting value (full vector of endogenous) +! shocks matrix of shocks (nexog x number of period) +! vcov covariance matrix of shocks (nexog x nexog) +! seed integer seed +! ysteady full vector of decision rule's steady +! dr structure containing matrices of derivatives (g_0, g_1,…) + +! output: +! res simulated results + +subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') + use iso_fortran_env + use iso_c_binding + use struct + use matlab_mex + use partitions + use simulation + 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) :: M_mx, options_mx, dr_mx, yhat_mx, epsilon_mx, udr_mx, tmp + type(pol), dimension(:), allocatable, target :: fdr, udr + integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nparticles, nys, nvar, nrestricted + real(real64), dimension(:), allocatable :: order_var, ys, ys_reordered, restrict_var_list, dyu + real(real64), dimension(:,:), allocatable :: yhat, e, ynext, ynext_all + type(pascal_triangle) :: p + type(uf_matching), dimension(:), allocatable :: matching + type(horner), dimension(:), allocatable :: h + integer :: i, d, j, m, n + character(kind=c_char, len=10) :: fieldname + + yhat_mx = prhs(1) + epsilon_mx = prhs(2) + dr_mx = prhs(3) + M_mx = prhs(4) + options_mx = prhs(5) + udr_mx = prhs(6) + + ! Checking the consistence and validity of input arguments + ! if (nrhs /= 5 .or. nlhs /= 1) then + ! call mexErrMsgTxt("Must have exactly 5 inputs and 1 output") + ! end if + if (nrhs /= 6 .or. nlhs /= 1) then + call mexErrMsgTxt("Must have exactly 5 inputs and 1 output") + end if + if (.not. (mxIsDouble(yhat_mx) .and. mxGetM(yhat_mx) >= 1 .and. mxGetN(yhat_mx) >= 1)) then + call mexErrMsgTxt("1st argument (yhat) should be a real vector") + end if + if (.not. (mxIsDouble(epsilon_mx) .and. mxGetM(epsilon_mx) >= 1 .or. mxGetN(epsilon_mx) == 1)) then + call mexErrMsgTxt("2nd argument (epsilon) should be a real vector") + end if + if (.not. mxIsStruct(dr_mx)) then + call mexErrMsgTxt("3rd argument (dr) should be a struct") + end if + if (.not. mxIsStruct(M_mx)) then + call mexErrMsgTxt("4th argument (M) should be a struct") + end if + if (.not. mxIsStruct(options_mx)) then + call mexErrMsgTxt("5th argument (options) should be a struct") + end if + if (.not. mxIsStruct(udr_mx)) then + call mexErrMsgTxt("6th argument (udr) should be a struct") + end if + + nstatic = get_int_field(M_mx, "nstatic") + npred = get_int_field(M_mx, "npred") + nboth = get_int_field(M_mx, "nboth") + nfwrd = get_int_field(M_mx, "nfwrd") + endo_nbr = nstatic + npred + nboth + nfwrd + exo_nbr = get_int_field(M_mx, "exo_nbr") + order = get_int_field(options_mx, "order") + nys = npred+nboth + nvar = nys + exo_nbr + + associate (order_var_mx => mxGetField(dr_mx, 1_mwIndex, "order_var")) + if (.not. (mxIsDouble(order_var_mx) .and. int(mxGetNumberOfElements(order_var_mx)) == endo_nbr)) then + call mexErrMsgTxt("Field dr.order_var should be a double precision vector with endo_nbr elements") + end if + allocate(order_var(endo_nbr)) + order_var = mxGetPr(order_var_mx) + end associate + + associate (ys_mx => mxGetField(dr_mx, 1_mwIndex, "ys")) + if (.not. (mxIsDouble(ys_mx) .and. int(mxGetNumberOfElements(ys_mx)) == endo_nbr)) then + call mexErrMsgTxt("Field dr.ys should be a double precision vector with endo_nbr elements") + end if + allocate(ys(endo_nbr), ys_reordered(endo_nbr)) + ys = mxGetPr(ys_mx) + ! Construct the reordered steady state + do i=1, endo_nbr + ys_reordered(i) = ys(int(order_var(i))) + end do + end associate + + associate (restrict_var_list_mx => mxGetField(dr_mx, 1_mwIndex, "restrict_var_list")) + if (.not. (mxIsDouble(restrict_var_list_mx))) then + call mexErrMsgTxt("Field dr.restrict_var_list should be a double precision vector") + end if + nrestricted = size(mxGetPr(restrict_var_list_mx)) + allocate(restrict_var_list(nrestricted)) + restrict_var_list = mxGetPr(restrict_var_list_mx) + end associate + + nparticles = int(mxGetN(yhat_mx)); + if (int(mxGetN(epsilon_mx)) /= nparticles) then + call mexErrMsgTxt("epsilon and yhat don't have the same number of columns") + end if + if (.not. (mxIsDouble(yhat_mx) .and. int(mxGetM(yhat_mx)) == npred + nboth)) then + call mexErrMsgTxt("yhat should be a double precision matrix with npred+nboth rows") + end if + if (.not. (mxIsDouble(epsilon_mx) .and. int(mxGetM(epsilon_mx)) == exo_nbr)) then + call mexErrMsgTxt("epsilon should be a double precision matrix with exo_nbr rows") + end if + + allocate(yhat(nys, nparticles), e(exo_nbr, nparticles), ynext(nrestricted, nparticles), ynext_all(endo_nbr, nparticles)) + yhat = reshape(mxGetPr(yhat_mx), [nys, nparticles]) + e = reshape(mxGetPr(epsilon_mx), [exo_nbr, nparticles]) + + allocate(h(0:order), fdr(0:order), udr(0:order)) + do i = 0, order + write (fieldname, '(a2, i1)') "g_", i + tmp = mxGetField(udr_mx, 1_mwIndex, trim(fieldname)) + if (.not. (c_associated(tmp) .and. mxIsDouble(tmp))) then + call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr") + end if + m = int(mxGetM(tmp)) + n = int(mxGetN(tmp)) + allocate(udr(i)%g(m,n), h(i)%c(endo_nbr, nvar**i)) + udr(i)%g(1:m,1:n) = reshape(mxGetPr(tmp), [m,n]) + end do + + ! Using the Horner algorithm to evaluate the decision rule at the chosen yhat and epsilon + allocate(dyu(nvar)) + do j=1,nparticles + dyu(1:nys) = yhat(:,j) + dyu(nys+1:) = e(:,j) + call eval(h, dyu, udr, endo_nbr, nvar, order) + ynext_all(:,j) = h(0)%c(:,1) + ys_reordered + do i=1,nrestricted + ynext(i,j) = ynext_all(int(restrict_var_list(i)),j) + end do + end do + + plhs(1) = mxCreateDoubleMatrix(int(size(restrict_var_list), mwSize), int(nparticles, mwSize), mxREAL) + mxGetPr(plhs(1)) = reshape(ynext, [size(ynext)]) + +end subroutine mexFunction \ No newline at end of file diff --git a/mex/sources/matlab_mex.F08 b/mex/sources/matlab_mex.F08 index 9c59871b9..865166317 100644 --- a/mex/sources/matlab_mex.F08 +++ b/mex/sources/matlab_mex.F08 @@ -250,6 +250,15 @@ module matlab_mat character(c_char), dimension(*), intent(in) :: fieldname end function mxGetField_internal + subroutine mxSetField_internal(pm, index, fieldname, pvalue) bind(c, name="mxSetField"//API_VER2) + use iso_c_binding + import :: mwIndex + type(c_ptr), intent(in), value :: pm + type(c_ptr), intent(in), value :: pvalue + integer(mwIndex), intent(in), value :: index + character(c_char), dimension(*), intent(in) :: fieldname + end subroutine mxSetField_internal + integer(c_int) function mxGetNumberOfFields(pm) bind(c, name="mxGetNumberOfFields"//API_VER) use iso_c_binding type(c_ptr), intent(in), value :: pm @@ -359,6 +368,14 @@ contains mxGetField = mxGetField_internal(pm, index-1, fieldname // c_null_char) end function mxGetField + subroutine mxSetField(pm, index, fieldname, pvalue) + type(c_ptr), intent(in) :: pm + type(c_ptr), intent(in) :: pvalue + integer(mwIndex), intent(in) :: index + character(kind=c_char, len=*), intent(in) :: fieldname + call mxSetField_internal(pm, index-1, fieldname // c_null_char, pvalue) + end subroutine mxSetField + type(c_ptr) function mxGetCell(pm, index) type(c_ptr), intent(in) :: pm integer(mwIndex), intent(in) :: index diff --git a/tests/particle/local_state_space_iteration_k_test.mod b/tests/particle/local_state_space_iteration_k_test.mod index bf09dc501..1a9c6c169 100644 --- a/tests/particle/local_state_space_iteration_k_test.mod +++ b/tests/particle/local_state_space_iteration_k_test.mod @@ -48,9 +48,15 @@ tStart1 = tic; for i=1:10000, ynext1 = local_state_space_iteration_2(yhat, epsil tStart2 = tic; for i=1:10000, ynext2 = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_); end, tElapsed2 = toc(tStart2); + +tStart3 = tic; [udr] = folded_to_unfolded_dr(dr, M_, options_); for i=1:10000, ynext3 = local_state_space_iteration_fortran(yhat, epsilon, dr, M_, options_, udr); end, tElapsed3 = toc(tStart3); + if max(max(abs(ynext1-ynext2))) > 1e-14 error('Inconsistency between local_state_space_iteration_2 and local_state_space_iteration_k') end +if max(max(abs(ynext1-ynext3))) > 1e-14 + error('Inconsistency between local_state_space_iteration_2 and local_state_space_iteration_fortran') +end if tElapsed1 1e-14 + error('Inconsistency between local_state_space_iteration_2 and local_state_space_iteration_fortran') +end + +if tElapsed32