function  [f,v] = dlqreig(a,b,q,r)
%
%DLQREIG  Discrete linear quadratic regulator solution
%  [f,v] = dlqreig(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
%
%  NOTE: For some problems it may be better to not balance the eigenvalue
%   problem.
%
%  ALGORITHM:  Eigenvector (with balancing)
%  WRITTEN BY:  Evan W. Anderson   December 1, 1994



%Calculate the eigenvector decomposition and find the indexes of the 
%stable and unstable eigenvalues

term = (b/r)*(b'/a');
[V,D] = eig([a+term*q   -term;  -(a'\q)     inv(a')]);
[E,I] = sort(abs(diag(D)));


%Now calculate the solution to the Riccati equation and the optimal
%decision rule

s = size(a,1);
v = V(s+1:2*s,I(1:s))/V(1:s,I(1:s));
f = (b'*v*b + r)\b'*v*a;