stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add iteration maximum to disclyap_fast.m prevent infinite loops

time-shift
Johannes Pfeifer 2015-07-22 15:08:25 +02:00
parent 5f7a8a6b0e
commit aee579ecc1
1 changed files with 21 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
matlab/partial_information/disclyap_fast.m
View File

@ -1,6 +1,14 @@
function X=disclyap_fast(G,V,tol,ch) function [X,exitflag]=disclyap_fast(G,V,tol,ch)
% function X=disclyap_fast(G,V,ch) % function X=disclyap_fast(G,V,ch)
% % Inputs:
% - G [double] first input matrix
% - V [double] second input matrix
% - tol [scalar] tolerance criterion
% - ch empty of non-empty if non-empty: check positive-definiteness
% Outputs:
% - X [double] solution matrix
% - exitflag [scalar] 0 if solution is found, 1 otherwise
%
% Solve the discrete Lyapunov Equation % Solve the discrete Lyapunov Equation
% X=G*X*G'+V % X=G*X*G'+V
% Using the Doubling Algorithm % Using the Doubling Algorithm
@ -11,7 +19,7 @@ function X=disclyap_fast(G,V,tol,ch)
% Joe Pearlman and Alejandro Justiniano % Joe Pearlman and Alejandro Justiniano
% 3/5/2005 % 3/5/2005
% Copyright (C) 2010-2012 Dynare Team % Copyright (C) 2010-2015 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -34,6 +42,7 @@ else
flag_ch = 1; flag_ch = 1;
end end
s=size(G,1); s=size(G,1);
exitflag=0;
%tol = 1e-16; %tol = 1e-16;
@ -41,13 +50,20 @@ P0=V;
A0=G; A0=G;
matd=1; matd=1;
while matd > tol iter=1;
while matd > tol && iter< 2000
P1=P0+A0*P0*A0'; P1=P0+A0*P0*A0';
A1=A0*A0; A1=A0*A0;
matd=max( max( abs( P1 - P0 ) ) ); matd=max( max( abs( P1 - P0 ) ) );
P0=P1; P0=P1;
A0=A1; A0=A1;
iter=iter+1;
end end
if iter==5000
X=NaN(P0);
exitflag=1;
return
end
clear A0 A1 P1; clear A0 A1 P1;
X=(P0+P0')/2; X=(P0+P0')/2;
@ -56,6 +72,7 @@ X=(P0+P0')/2;
if flag_ch==1 if flag_ch==1
[C,p]=chol(X); [C,p]=chol(X);
if p ~= 0 if p ~= 0
exitflag=1;
error('X is not positive definite') error('X is not positive definite')
end end
end end