2005-02-18 20:54:39 +01:00
|
|
|
function AA = generalized_cholesky(A);
|
2008-09-26 18:55:26 +02:00
|
|
|
%function AA = generalized_cholesky(A);
|
|
|
|
%
|
|
|
|
% Calculates the Gill-Murray generalized choleski decomposition
|
|
|
|
% Input matrix A must be non-singular and symmetric
|
|
|
|
|
|
|
|
% Copyright (C) 2003 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 <http://www.gnu.org/licenses/>.
|
2005-02-18 20:54:39 +01:00
|
|
|
|
|
|
|
n = rows(A);
|
|
|
|
R = eye(n);
|
|
|
|
E = zeros(n,n);
|
|
|
|
norm_A = max(transpose(sum(abs(A))));
|
|
|
|
gamm = max(abs(diag(A)));
|
|
|
|
delta = max([eps*norm_A;eps]);
|
|
|
|
|
2008-10-04 23:22:08 +02:00
|
|
|
for j = 1:n
|
2005-02-18 20:54:39 +01:00
|
|
|
theta_j = 0;
|
2008-10-04 23:22:08 +02:00
|
|
|
for i=1:n
|
2005-02-18 20:54:39 +01:00
|
|
|
somme = 0;
|
2010-01-05 11:46:10 +01:00
|
|
|
for k=1:i-1
|
2008-10-04 23:22:08 +02:00
|
|
|
somme = somme + R(k,i)*R(k,j);
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
R(i,j) = (A(i,j) - somme)/R(i,i);
|
2008-10-04 23:22:08 +02:00
|
|
|
if (A(i,j) -somme) > theta_j
|
|
|
|
theta_j = A(i,j) - somme;
|
|
|
|
end
|
|
|
|
if i > j
|
|
|
|
R(i,j) = 0;
|
|
|
|
end
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
somme = 0;
|
2010-01-05 11:46:10 +01:00
|
|
|
for k=1:j-1
|
2005-02-18 20:54:39 +01:00
|
|
|
somme = somme + R(k,j)^2;
|
2008-10-04 23:22:08 +02:00
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
phi_j = A(j,j) - somme;
|
2008-10-04 23:22:08 +02:00
|
|
|
if j+1 <= n
|
2005-02-18 20:54:39 +01:00
|
|
|
xi_j = max(abs(A((j+1):n,j)));
|
2008-10-04 23:22:08 +02:00
|
|
|
else
|
2005-02-18 20:54:39 +01:00
|
|
|
xi_j = abs(A(n,j));
|
2008-10-04 23:22:08 +02:00
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
beta_j = sqrt(max([gamm ; (xi_j/n) ; eps]));
|
2008-10-04 23:22:08 +02:00
|
|
|
if all(delta >= [abs(phi_j);((theta_j^2)/(beta_j^2))])
|
2005-02-18 20:54:39 +01:00
|
|
|
E(j,j) = delta - phi_j;
|
2008-10-04 23:22:08 +02:00
|
|
|
elseif all(abs(phi_j) >= [((delta^2)/(beta_j^2));delta])
|
2005-02-18 20:54:39 +01:00
|
|
|
E(j,j) = abs(phi_j) - phi_j;
|
2008-10-04 23:22:08 +02:00
|
|
|
elseif all(((theta_j^2)/(beta_j^2)) >= [delta;abs(phi_j)])
|
2005-02-18 20:54:39 +01:00
|
|
|
E(j,j) = ((theta_j^2)/(beta_j^2)) - phi_j;
|
2008-10-04 23:22:08 +02:00
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
R(j,j) = sqrt(A(j,j) - somme + E(j,j));
|
2008-10-04 23:22:08 +02:00
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
AA = transpose(R)*R;
|