function X=disclyap_fast(G,V,ch) Solve the discrete Lyapunov Equation X=G*X*G'+V Using the Doubling Algorithm If ch is defined then the code will check if the resulting X is positive definite and generate an error message if it is not Joe Pearlman and Alejandro Justiniano 3/5/2005
0001 function X=disclyap_fast(G,V,tol,ch) 0002 % function X=disclyap_fast(G,V,ch) 0003 % 0004 % Solve the discrete Lyapunov Equation 0005 % X=G*X*G'+V 0006 % Using the Doubling Algorithm 0007 % 0008 % If ch is defined then the code will check if the resulting X 0009 % is positive definite and generate an error message if it is not 0010 % 0011 % Joe Pearlman and Alejandro Justiniano 0012 % 3/5/2005 0013 0014 % Copyright (C) 2010-2011 Dynare Team 0015 % 0016 % This file is part of Dynare. 0017 % 0018 % Dynare is free software: you can redistribute it and/or modify 0019 % it under the terms of the GNU General Public License as published by 0020 % the Free Software Foundation, either version 3 of the License, or 0021 % (at your option) any later version. 0022 % 0023 % Dynare is distributed in the hope that it will be useful, 0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0026 % GNU General Public License for more details. 0027 % 0028 % You should have received a copy of the GNU General Public License 0029 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0030 0031 if nargin <= 3 || isempty( ch ) == 1 0032 flag_ch = 0; 0033 else 0034 flag_ch = 1; 0035 end 0036 s=size(G,1); 0037 0038 %tol = 1e-16; 0039 0040 P0=V; 0041 A0=G; 0042 0043 matd=1; 0044 while matd > tol 0045 P1=P0+A0*P0*A0'; 0046 A1=A0*A0; 0047 matd=max( max( abs( P1 - P0 ) ) ); 0048 P0=P1; 0049 A0=A1; 0050 end 0051 clear A0 A1 P1; 0052 0053 X=(P0+P0')/2; 0054 0055 % Check that X is positive definite 0056 if flag_ch==1 0057 [C,p]=chol(X); 0058 if p ~= 0 0059 error('X is not positive definite') 0060 end 0061 end