From f392c78644c0d1d29a0bd6215686b47526934ab5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 21 Dec 2023 16:32:22 +0100 Subject: [PATCH] Drop unused riccati_update MEX --- meson.build | 5 -- mex/sources/riccati_update/mexFunction.f08 | 98 ---------------------- tests/riccatiupdate.m | 87 ------------------- 3 files changed, 190 deletions(-) delete mode 100644 mex/sources/riccati_update/mexFunction.f08 delete mode 100644 tests/riccatiupdate.m diff --git a/meson.build b/meson.build index 6855c0c08..b0f64c019 100644 --- a/meson.build +++ b/meson.build @@ -346,10 +346,6 @@ shared_module('logarithmic_reduction', [ 'mex/sources/logarithmic_reduction/mexF shared_module('disclyap_fast', [ 'mex/sources/disclyap_fast/disclyap_fast.f08' ] + mex_blas_fortran_iface, kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ]) -# TODO: Same remark as A_times_B_kronecker_C -shared_module('riccati_update', [ 'mex/sources/riccati_update/mexFunction.f08' ] + mex_blas_fortran_iface, - kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ]) - qmc_sequence_src = [ 'mex/sources/sobol/qmc_sequence.cc', 'mex/sources/sobol/sobol.f08' ] # Hack for statically linking libgfortran @@ -1811,7 +1807,6 @@ mod_and_m_tests = [ 'solver-test-functions/wood.m',] }, { 'test' : [ 'cyclereduction.m' ] }, { 'test' : [ 'logarithmicreduction.m' ] }, - { 'test' : [ 'riccatiupdate.m' ] }, { 'test' : [ 'kalman/likelihood/test_kalman_mex.m' ] }, { 'test' : [ 'contribs.m' ], 'extra' : [ 'sandbox.mod', diff --git a/mex/sources/riccati_update/mexFunction.f08 b/mex/sources/riccati_update/mexFunction.f08 deleted file mode 100644 index 7443a604c..000000000 --- a/mex/sources/riccati_update/mexFunction.f08 +++ /dev/null @@ -1,98 +0,0 @@ -! Copyright © 2022-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 . - -! Implements Ptmp = T*(P-K*Z*P)*transpose(T)+Q where -! P is the (r x r) variance-covariance matrix of the state vector -! T is the (r x r) transition matrix of the state vector -! K is the (r x n) gain matrix -! Z is the (n x r) matrix linking observable variables to state variables -! Q is the (r x r) variance-covariance matrix of innovations in the state equation -! and accounting for different properties: -! P is a (symmetric) positive semi-definite matrix -! T can be triangular - -subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') - use matlab_mex - use blas - implicit none (type, external) - - type(c_ptr), dimension(*), intent(in), target :: prhs - type(c_ptr), dimension(*), intent(out) :: plhs - integer(c_int), intent(in), value :: nlhs, nrhs - - real(real64), dimension(:,:), pointer, contiguous :: P, T, K, Z, Q, Pnew - real(real64), dimension(:,:), allocatable :: tmp1, tmp2 - integer :: i, n, r - character(kind=c_char, len=2) :: num2str - - ! 0. Checking the consistency and validity of input arguments - if (nrhs /= 5_c_int) then - call mexErrMsgTxt("Must have 5 input arguments") - end if - if (nlhs > 1_c_int) then - call mexErrMsgTxt("Too many output arguments") - end if - - do i=1,5 - if (.not. (c_associated(prhs(i)) .and. mxIsDouble(prhs(i)) .and. & - (.not. mxIsComplex(prhs(i))) .and. (.not. mxIsSparse(prhs(i))))) then - write (num2str,"(i2)") i - call mexErrMsgTxt("Argument " // trim(num2str) // " should be a real dense matrix") - end if - end do - - r = int(mxGetM(prhs(1))) ! Number of states - n = int(mxGetN(prhs(3))) ! Number of observables - - if ((r /= mxGetN(prhs(1))) & ! Number of columns of P - &.or. (r /= mxGetM(prhs(2))) & ! Number of lines of T - &.or. (r /= mxGetN(prhs(2))) & ! Number of columns of T - &.or. (r /= mxGetM(prhs(3))) & ! Number of lines of K - &.or. (n /= mxGetM(prhs(4))) & ! Number of lines of Z - &.or. (r /= mxGetN(prhs(4))) & ! Number of columns of Z - &.or. (r /= mxGetM(prhs(5))) & ! Number of lines of Q - &.or. (r /= mxGetN(prhs(5))) & ! Number of columns of Q - ) then - call mexErrMsgTxt("Input dimension mismatch") - end if - - ! 1. Storing the relevant information in Fortran format - P(1:r,1:r) => mxGetPr(prhs(1)) - T(1:r,1:r) => mxGetPr(prhs(2)) - K(1:r,1:n) => mxGetPr(prhs(3)) - Z(1:n,1:r) => mxGetPr(prhs(4)) - Q(1:r,1:r) => mxGetPr(prhs(5)) - - plhs(1) = mxCreateDoubleMatrix(int(r, mwSize), int(r, mwSize), mxREAL) - Pnew(1:r, 1:r) => mxGetPr(plhs(1)) - - ! 2. Computing the Riccati update of the P matrix - allocate(tmp1(r,r), tmp2(r,r)) - ! Pnew <- Q - Pnew = Q - ! tmp1 <- K*Z - call matmul_add("N", "N", 1._real64, K, Z, 0._real64, tmp1) - ! tmp2 <- P - tmp2 = P - ! tmp2 <- tmp2 - tmp1*P - call matmul_add("N", "N", -1._real64, tmp1, P, 1._real64, tmp2) - ! tmp1 <- T*tmp2 - call matmul_add("N", "N", 1._real64, T, tmp2, 0._real64, tmp1) - ! Pnew <- tmp1*T' + Pnew - call matmul_add("N", "T", 1._real64, tmp1, T, 1._real64, Pnew) - -end subroutine mexFunction diff --git a/tests/riccatiupdate.m b/tests/riccatiupdate.m deleted file mode 100644 index f293c2fe0..000000000 --- a/tests/riccatiupdate.m +++ /dev/null @@ -1,87 +0,0 @@ -source_dir = getenv('source_root'); -addpath([source_dir filesep 'matlab']); - -dynare_config; - -testFailed = 0; - -skipline() -disp('*** TESTING: riccatiupdate.m ***'); - -t0 = clock; - -% Set the number of experiments for time measurement -N = 5000; -% Set the dimension of the problem to be solved. -r = 50; -n = 100; -tol = 1e-15; -% Set the input arguments -% P, Q: use the fact that for any real matrix A, A'*A is positive semidefinite -P = rand(n,r); -P = P'*P; -Q = rand(n,r); -Q = Q'*Q; -K = rand(r,n); -Z = rand(n,r); -T = rand(r,r); -% Computing an upperbound for the norm the updated variance-covariance matrix -ub = norm(T,1)^2*norm(P,1)*(1+norm(K*Z,1))+norm(Q,1); -% Weighting the P and Q matrices to keep the norm of the variance-covariance matrix below 1 -P = 0.5*P/ub; -Q = 0.5*Q/ub; - -% 1. Update the state vairance-covariance matrix with Matlab -tElapsed1 = 0.; -tic; -for i=1:N - Ptmp_matlab = T*(P-K*Z*P)*transpose(T)+Q; -end -tElapsed1 = toc; -disp(['Elapsed time for the Matlab Riccati update is: ' num2str(tElapsed1) ' (N=' int2str(N) ').']) - -% 2. Update the state varance-covariance matrix with the mex routine -tElapsed2 = 0.; -Ptmp_fortran = P; -try - tic; - for i=1:N - Ptmp_fortran = riccati_update(P, T, K, Z, Q); - end - tElapsed2 = toc; - disp(['Elapsed time for the Fortran Riccati update is: ' num2str(tElapsed2) ' (N=' int2str(N) ').']) - R = norm(Ptmp_fortran-Ptmp_matlab,1); - if (R > tol) - testFailed = testFailed+1; - dprintf('The Fortran Riccati update is wrong') - end -catch - testFailed = testFailed+1; - dprintf('Fortran Riccati update failed') -end - -% Compare the Fortran and Matlab execution time -if tElapsed1 0)