function [P1,status,NumIter,TestValue] = RiccatiSolver(F,H,Q,R,P1,MaxIter,ConvValue)
% RiccatiSolver: Brute force solver of the discrete algebraic Riccati equation when Q and R
%                can both be singular.
%
% USAGE:
%
%       [P1,status,NumIter,TestValue] = RiccatiSolver(F,H,Q,R,P1,MaxIter,ConvValue)
%
% where the discrete algebraic Riccati equation is expressed as:
%
%         P1 = F*P1*F' - F*P1*H*inverse( H'*P1*H + R )*H'*P1*F' + Q (1)
%
% The RiccatiSolver function iterates on until the largest absolute change in any element of P1 from iteration i to i+1
% is smaller than or equal to the tolerance value or the dare function has converged to a unique solution. For each new
% P1 matrix it uses Schur decompositions of the matrices P1 into N*Lambda*N' and  B'*Lambda*B+R, where B = N'*H. These
% decompositions make it possible to use the dare algorithm for Lambda.
%
% REQUIRED INPUT:  F (rxr) transition matrix from the state equation.
%
%                  H rxn) matrix where n is the number of observed variables and r the number of
%                         state variables. This matrix is part of the measurement equation of the
%                         state-space representation for the Kalman filter.
%
%                  Q (rxr) covariance matrix for the shocks to the state equation.
%
%                  R (nxn) covariance matrix for the measurement errors in the measurement equation.
%
%                  P1 (rxr) matrix with initial values for the P1 matrix.
%
%                  MaxIter (integer) with the maximum number of iterations.
%
%                  ConvValue (scalar), a positive real number measuring the tolerance value.
%
% REQUIRED OUTPUT: P1 (rxr) matrix with the numerical solution to the Riccati equation.
%
%                  status (integer) takes the value 0 if the solution method converged and 1 otherwise. It takes
%                         on the value 2 in the special case that infinite or NaN's are found in P1 during the
%                         iterations.
%
%                  NumIter (integer) with the number of iterations used to obtain convergence.
%
% OPTIONAL OUTPUT: TestValue (scalar), with value of the convergence criterion, the matlab norm function
%                         of the difference between the last two P1 matrices or (if dare converged) the
%                         absolute value of the relative residual from dare.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: October 24, 2006.
%                        This version: December 5, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA 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.
%
%      This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 26-10-2006: Added the drawnow call in the loop.
%
% * 31-10-2006: Updated the documentation.
%
% * 8-11-2006:  Added the dare step of the algorithm and updated the documentation.
%
% * 9-11-2006:  Took the name change of dare to YADAdare into account.
%
% * 20-11-2006: Added the test for infinite or NaN's at the start of each iteration.
%
% * 23-11-2006: Added additional tests for Inf and NaN's before calling the SchurDecomposition function.
%
% * 2-11-2007:  Updated the documentation.
%
% * 23-5-2008:  Updated the documentation.
%
% * 21-7-2008:  The function now checks if it should use the internal dare function or the external.
%
% * 5-12-2008:  Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

status = 1;
NumIter = MaxIter;
n = size(H,2);
%
if exist('YADAdare')==2;
   UseYADAdare = 1;
else;
   UseYADAdare = 0;
end;
%
for i=1:MaxIter;
   if (sum(sum(isinf(P1)))>0)|(sum(sum(isnan(P1)))>0);
      if nargout==4;
         TestValue = 0;
      end;
      NumIter = i;
      status = 2;
      break;
   end;
   P = P1;
   if rank((H'*P*H)+R)<n;
      S = SchurDecomposition((H'*P*H)+R);
   else;
      S = eye(n);
   end;
   P1 = (F*P*F')-(F*P*H*S*inv(S'*((H'*P*H)+R)*S)*S'*H'*P*F')+Q;
   %
   % try the new P1 matrix for NaN's and Inf
   %
   if (sum(sum(isinf(P1)))>0)|(sum(sum(isnan(P1)))>0);
      if nargout==4;
         TestValue = 0;
      end;
      NumIter = i;
      status = 2;
      break;
   end;
   %
   % apply a Schur decomposition of P1
   %
   [N,Lambda] = SchurDecomposition(P1);
   %
   % try to solve the equation taking reduced rank into account
   %
   A = N'*F'*N;
   B = N'*H;
   C = N'*Q*N;
   Rn = R;
   %
   % try the (B'*Lambda*B)+Rn matrix for Inf and NaN's
   %
   if (sum(sum(isinf((B'*Lambda*B)+Rn)))>0)|(sum(sum(isnan((B'*Lambda*B)+Rn)))>0);
      if nargout==4;
         TestValue = 0;
      end;
      NumIter = i;
      status = 2;
      break;
   end;
   %
   % check the rank of B'*Lambda*B+Rn
   %
   if rank((B'*Lambda*B)+Rn)<size((B'*Lambda*B)+Rn,1);
      D = SchurDecomposition((B'*Lambda*B)+Rn);
   else;
      D = eye(size((B'*Lambda*B)+Rn,1));
   end;
   if UseYADAdare==1;
      [X,L,G,rep] = YADAdare(A,B*D,C,D'*Rn*D,'report');
   else;
      [X,L,G,rep] = dare(A,B*D,C,D'*Rn*D,'report');
   end;
   %
   % check if we have convergence
   %
   if (norm(P1-P)<=ConvValue)|((isnan(rep)==0)&(rep>-1));
      if (isnan(rep)==0)&(rep>-1);
         %
         % the dare algorithm converged
         %
         P1 = N*X*N';
         if rank((H'*P1*H)+R)<n;
            S = SchurDecomposition((H'*P1*H)+R);
         else;
            S = eye(n);
         end;
         P2 = (F*P1*F')-(F*P1*H*S*inv(S'*((H'*P1*H)+R)*S)*S'*H'*P1*F')+Q;
      end;
      status = 0;
      NumIter = i;
      break;
   end;
   %
   % update the screen to make sure it doesn't look like it has frozen
   %
   drawnow;
end;

if (nargout==4)&(status~=2);
   TestValue = norm(P1-P);
   if (isnan(rep)==0)&(rep>-1);
      TestValue = abs(rep);
      if (sum(sum(isinf(P1)))==0)&(sum(sum(isnan(P1)))==0)&(sum(sum(isinf(P2)))==0)&(sum(sum(isnan(P2)))==0);
         TestValue = norm(P1-P2);
      end;
   end;
end;

%
% end of RiccatiSolver.m
%
