function  [f,v] = dlqrmsign(a,b,q,r,tol,maxit)
%
%DLQRMSIGN  Discrete linear quadratic regulator solution
%  [f,v] = dlqrmsign(a,b,q,r) solves the Riccati equation
%  associated with the real discrete time optimal linear regulator
%  problem of the form
%
%        min     sum {x'qx + u'ru}
%        st      x(t+1) = ax(t) + bu(t)
%
%  Returned are the matrices f and v such that the optimal decision rule
%  is -fx and the value function is x'vx  

%   tol and maxit are optional
%  parameters.  If they both are not inputted then they are given 
%  default values.  See the m-file msign.m for a description of
%  tol and maxit. 

%
%  ALGORITHM: Matrix sign - symmetry is exploited and the det weighting
%     scheme is used
%  REFERENCES:
%     Gardiner and Laub(1986)
%     Laub(1991)
%  WRITTEN BY:  Evan W. Anderson       Dec  1, 1994
%  DOCUMENTATION LAST MODIFIED: EWA    Dec 29, 1998


s = size(a,1);
if nargin < 6, 
  tol = 1e-15; 
  maxit = 1e4; 
end 

%Calculate the Matrix Sign  J*MBAR

N = [a zeros(s); -q  eye(s)];
L = [eye(s) b*(r\b'); zeros(s)  a'];
MBAR = (N-L)\(N+L);
[Y, INFO]= msign([MBAR(s+1:2*s,:); -MBAR(1:s,:)],tol,maxit);


%Calculate the solution to the Riccati equation and optimal decision rule
%Let Z = J*Y . Then the solution is:
%
%            (  Z12   )   \     (Z11 + I)
%     v = -  ( Z22 + I)    \    (Z21)

v = [Y(s+1:2*s,s+1:2*s); -Y(1:s,s+1:2*s) - eye(s)]\  ...
         [-Y(s+1:2*s,1:s)+eye(s); Y(1:s,1:s) ]; 
f = (b'*v*b + r)\b'*v*a;