Drop unused riccati_update MEX
Sébastien Villemot
parent
fbb89dc190
commit
f392c78644
|
|
@ -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',
|
||||
|
|
|
|||
|
|
@ -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 <https://www.gnu.org/licenses/>.
|
||||
|
||||
! 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
|
||||
|
|
@ -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<tElapsed2
|
||||
skipline()
|
||||
dprintf('Matlab Riccati update is %5.2f times faster than its Fortran counterpart.', tElapsed2/tElapsed1)
|
||||
skipline()
|
||||
else
|
||||
skipline()
|
||||
dprintf('Fortran Riccati update is %5.2f times faster than its Matlab counterpart.', tElapsed1/tElapsed2)
|
||||
skipline()
|
||||
end
|
||||
|
||||
% Compare results after multiple calls
|
||||
N = 50;
|
||||
disp(['After 1 update using the Riccati formula, the norm-1 discrepancy is ' num2str(norm(Ptmp_fortran-Ptmp_matlab,1)) '.']);
|
||||
for i=2:N
|
||||
Ptmp_matlab = T*(Ptmp_matlab-K*Z*Ptmp_matlab)*transpose(T)+Q;
|
||||
Ptmp_fortran = riccati_update(Ptmp_fortran, T, K, Z, Q);
|
||||
disp(['After ' int2str(i) ' updates using the Riccati formula, the norm-1 discrepancy is ' num2str(norm(Ptmp_fortran-Ptmp_matlab,1)) '.'])
|
||||
end
|
||||
|
||||
t1 = clock;
|
||||
|
||||
fprintf('\n*** Elapsed time (in seconds): %.1f\n\n', etime(t1, t0));
|
||||
|
||||
quit(testFailed > 0)
|
||||
Loading…
Reference in New Issue