function [F,P]=ricroots(beta,A,B,Q,R,W);
%RICROOTS uses the Newton-Raphson routine to find the fixed point of the 
%  Riccati system:
%
%        F=(R+beta*B'*P*B)\(beta*B'*P*A+W')
%        P=Q+beta*A'*P*A-(beta*A'*P*B+W)*F
% 
%  the decision rule and Riccati equation for the optimal linear regulator
%  problem:
%          sum {beta^t [x'Qx + u'Ru +2x'Wu] }
%       
%  subject to 
%		x[t+1] = Ax[t] + Bu[t] 
%
%  where x is the mx1 vector of states, u is the kx1 vector of controls,
%  A is mxm, B is mxk, Q is mxm, R is kxk, W is mxk.
%                
%  INPUTS: beta, A, B, Q, R, W  (or just beta, A, B, Q, R if there are no
%          cross products, W=0).  Set beta=1, if there is no discounting.
%
%  OUTPUTS: 
%    (i) F= (R+beta*B'*p*B)\(beta*B'*p*A+W'); the decision rule: u=-F x
%   (ii) P= beta*A'*P*A + Q -(beta*A'*P*B+W)*f, the solution to the matrix
%           Riccati equation
%
%  NOTES: see also DOUBLE, DRICC and OLRP.
%
[n,m]=size(B);
s=(m+n)*n;
P=-eye(n);
F0= (R+beta*B'*P*B)\(beta*B'*P*A+W');
dd=1;
it=1;
while (dd>1e-10 & it<=5);
  P=  beta*A'*P*A + Q -(beta*A'*P*B+W)*F0;
  F1= (R+beta*B'*P*B)\(beta*B'*P*A+W');
  dd=max(max(abs( (F1-F0)./max(ones(abs(sign(F0))),abs(F0)) )));
  it=it+1;
  F0=F1;
end;
F=F1;
x0=[F(:);P(:)];
for i=1:100;
  tf=(R+beta*B'*P*B)*F-beta*B'*P*A-W';
  tp=P-Q-beta*A'*P*A+(beta*A'*P*B+W)*F;
  h=[tf(:);
     tp(:)];
  hj=[kron(eye(n),R+beta*B'*P*B)  kron(beta*(F'*B'-A'),B');
      kron(eye(n),beta*A'*P*B+W)  eye(n^2)+kron(beta*(F'*B'-A'),A')];
  x=x0-hj\h;
  if max(abs( (x-x0)./max(ones(abs(sign(x))),abs(x0)) ))<1e-10; break; end;
  x0=x;
  F(:)=x0(1:n*m);
  P(:)=x0(n*m+1:s);
end;
if i>=100; disp('Iteration limit of 100 reached in RICROOTS'); end;
